Copula Graphical Models for Wind Resource Estimation

1 Introduction

This paper addresses wind resource assessment: the problem of determining if there will be enough wind in the ideal speed range that will endure at a potential wind farm or "site", over a 20+ year timespan. The end-to-end pipeline of a resource assessment service spans from automatic site-neighbor data extraction from public, online sources (ASOS database), through site-neighbor data synchronization in preparation for generative modeling, modeling, backcast (where historical data at neighboring sites is passed through a model to obtain predictions at the site) to estimation of the industry standard Weibull distribution from the derived predictions. Our focus in this paper is on the single most critical factor in assessment: achieving the most accurate backcast while incurring minimal financial expense. This implies integrating geographically proximal public wind data sources and building models that better represent the data (accuracy) while concurrently reducing the duration of anemometer sensing during the assessment period (expense).

The primary contribution of this paper is the use of multivariate copulas for modeling multiple joint distributions of wind speeds at the site and a publicly available neighboring wind source. We construct a n-dimensional Gaussian copula and multiple copula graphical models to enhance the quality of the prediction site distribution. This modeling step is embedded within a widespread methodology called Measure-Corelate-Predict (MCP) [Gross and Phelan, 2006; Bass et al., 2000; Bailey et al., 1997; Lackner et al., 2008].

For demonstration we use speed and direction data from an actual site in the state of Massachusetts where we have assembled data from anemometer sensing carried over for a period of two years. We compare our models with multiple regression methods, where it achieves higher accuracy with Una-May O'Reilly CSAIL, MIT Cambridge, MA unamay@csail.mit.edu less sensing data - sometimes with only 3 months. The industry standard method, multiple regression, achieves a reasonable accuracy with 8 months of data, an industry standard period. Thus we achieve better accuracy at a lower cost.

We proceed by describing MCP while introducing notation in Section 2. Section 3 describes the real wind resource estimation scenario and the dataset we utilized throughout this paper to demonstrate our methods. Section 4 describes the copula modeling. Section 5 is the demonstration. We intentionally reference related work throughout the paper, in context of discussion.

2 Measure-Corelate-Predict (MCP )

We consider wind resource estimation derived by a methodology known as Measure-Correlate-Predict or MCP. In terms of notation, the wind at a particular location is characterized by speed denoted by x and direction θ. Wind speed is measured by anemometers and wind direction is measured by wind vanes. The 360o direction is split into multiple bins with a lower limit (θl) and upper limit (θu). We give an index value of J = 1 . . . j for the directional bin. We represent the wind speed measurement at the test site (where wind resource needs to be estimated) with y and the other sites (for whom the long term wind resource is available) as x and index these other sites with M = 1 . . . m.

The three steps of MCP are:

MEASURE Short term sensing measurements on the site are collected. This is denoted by Y = {ytk . . . ytn }. Measurements can be collected using anemometers on the site, a newly-constructed meteorological tower, or even remote sensing technologies such as sonar or lidar. Different measurement techniques incur different costs that dictate their feasibility for different projects. Measurements from nearby sites for the same period are gathered. These sites, called historical sites, have additional data for the past 10-20 years. These are denoted by X = {x1...m tk...tn } where each xitk...tn corresponds to data from one historical site and m denotes the total number of historical sites. Historical data that is not simultaneous in time to the site observations used in modeling will be used in the PREDICT step.

CORRELATE A single directional model is first built correlating the wind directions observed at the site with simultaneous historical site wind directions. Next, for each directional interval, called a (directional) bin, of a 360◦ radius, a model is built correlating the wind speeds at the site with simultaneous speeds at the historical sites, i.e. Yti = fθj (x1...m ti ) where k ≤ i ≤ n. The data available from the site at this stage is expected to be sparse and noisy.

PREDICT To obtain an accurate estimation of long term wind conditions at the site, we first divide the data from the historic sites (which is not simultaneous in time to the site observations used in modeling) into subsets that correspond to a directional bin. Prediction of the long term site conditions follows two steps:

A : We use the model we developed for that direction fθj and the data from the historic sites corresponding to this direction x1...m t1...tk-1 |θj to predict what the wind speed Yp = yt1...tk-1 at the site would be. For a new observation x we have to predict y. For this we form the conditional first by

Our predicted ˆ y maximizes this conditional probability

Note that the term in the denominator of eq.( 1) remains constant, hence for the purposes of finding the optimum we can ignore its evaluation. We simply evaluate this conditional for the entire range of Y in discrete steps and pick the value of y ∈ Y that maximizes the conditional.

B : With the predictions Yp, from A above, we estimate parameters for a Weibull distribution. This distribution is our answer to the wind resource assessment problem. We generate a distribution for each directional bin.

The goal is to generate a predicted long term wind speed distribution in each direction which will be as close as possible to the real (as yet unexperienced) distribution. The result fromMCP, i.e. the statistical distribution in each bin, is then used to estimate the energy which can be expected from a wind turbine, given the power curve supplied by its manufacturer. This calculation can be extended over an entire farm if wake interactions among the turbines are taken into account. See [Wagner et al., 2011] for more details. Note that distribution not only captures the mean, but also variance in this speed. This is critical for assessment of long term wind resource and the long term energy estimate.

A variety of methods are developed in [Rogers et al., 2005] to evaluate the accuracy of the predicted wind speed distribution. One method measures the accuracy in terms of ratios between true and actual parameters of the Weibull distribution. That is, true shape versus estimated shape and true scale versus estimated scale. To completely capture any possible inaccuracy in the predicted distribution, we measure a symmetric Kullback-Leibler distance. It is important to note that this measure is different than the mean-squared error or mean-absolute error which measure the accuracy in terms of difference between each predicted value and the true observation. Methods that minimize these errors would not necessarily accurately express how close the approximation is to the true distribution. As a measure of predictive accuracy we compare the final estimated Weibull distribution to the ground truth distribution using Kullback-Leibler (KL) divergence. The lower this value, the more accurate the prediction:

KL divergence derives the distance between two probability distributions:

For baseline comparison, we also developed a linear regression model which is used quite extensively in wind resource assessment [Bass et al., 2000; Rogers et al., 2005].

We now proceed to describe the our machine learning approaches for wind resource assessment.

3 A real world scenario and dataset

To evaluate and compare our different algorithms, we acquired wind data collected using anemometers from the rooftop of Museum of Science in Boston where a wind vane is also installed. These anemometers are inexpensive and consequently noisy. The museum is located amongst buildings, a river and is close to a harbor as shown in Figure 1. This provides us with a site that is topographically challenging. At this location we have approximately 2 years worth of data collected at a frequency of 1 sample/second with 10 minute averages stored in a separate database. To derive the wind resource assessment we train using data from the first year. This data is split into three datasets we call D3, D6 and D8. The split D3 has data for 3 months. The split D6 has 3 additional months for a total of 6 and D8 has yet 2 more months for a total of 8. We divide each dataset and the second year's dataset (to serve as our test data, i.e. ground truth) further into 12 directional bins of equal sizes starting at compass point North (0◦ ).

We use airport wind data from the public ASOS (Automated Surface Observing System) database for sources of neighboring-site data. This data is regularly accessed by the wind industry for correlation purposes. The airports' locations are shown in Figure 1 (right).

4 Multivariate copulas

Previous modeling techniques assume a Gaussian distribution for wind speed and direction variables and a Gaussian joint distribution. It is arguable however that Gaussian distributions do not accurately represent the wind speed distributions. In fact, conventionally a univariate Weibull distribution [Burton et al., 2001] is used to parametrically describe wind sensor measurements. A Weibull distribution is likely also chosen for its flexibility because it can express any one of multiple distributions, including Rayleigh or Gaussian.

To the best of our knowledge, however, joint density functions for non-Gaussian distributions have not been estimated for wind resource assessment. In this paper, to build a multivariate model from marginal distributions which are not all Gaussian, we exploit copula functions. A copula framework provides a means of inference after modeling a multivariate joint distribution from training data.

Because copula estimation is less well known, we now briefly review copula theory. We will then describe how we construct the individual parametric distributions which are components of a copula and then, how we couple them to form a multivariate density function. Finally, we present our approach to predict the value of y given x1...m.

A copula function C(u1, . . . um+1; θ) with parameter θ represents a joint distribution function for multiple uniform random variables U1 . . . Um+1 such that

Let U1 . . .Um represent the cumulative distribution functions (CDF) for variables x1, . . . xm and Um+1 represent the CDF for y. Hence the copula represents the joint distribution function of C(F(x1) . . . F(xm), F(y)), where Ui = F(xi). According to Sklar's theorem any copula function taking marginal distributions F(xi) as its arguments, defines a valid joint distribution with marginals F(xi). Thus we are able to construct the joint distribution function for x1 . . . xm, y:

The joint probability density function (PDF) is obtained by taking the m+ 1th order derivative of eqn. (6)

where c(.) is the copula density. Thus the joint density function is a weighted version of independent density functions, where the weight is derived via copula density.

4.1 Gaussian copula

First we consider a multivariate Gaussian copula to form a statistical model for our variables given by

where FG is the CDF of multivariate normal with zero mean vector and Σ as covariance and F -1 is the inverse of the standard normal.

Estimation of parameters: There are two sets of parameters to estimate. The first set of parameters for the multivariate Gaussian copula is Σ. The second set, denoted by Ψ = {ψ, ψy} are the parameters for the marginals of x, y. Given N i.i.d observations of the variables x, y, the loglikelihood function is:

Parameters Ψ are estimated via[Iyengar, 2011]

A variety of algorithms are available in literature to estimate the MLE in eq. (10). We refer users to [Iyengar, 2011] for a thorough discussion of estimation methods. For more details about the copula theory readers are referred to [Nelsen, 2006].

4.2 Vine models

Before giving details of how to construct a vine, we present three examples of how to derive the factorization of a multivariate probability distribution in terms of bivariate copulas. Consider first only two variables x1 and x2. The joint density, f12(x1, x2), can be factorized in two ways. First, using the chain rule:

• (a) f(x1, . . . , x5) = (c1,2c2,3c3,4c4,5)(c1,3|2c2,4|3c3,5|4)(c1,4|2,3c2,5|3,4)(c1,5|2,3,4) ∏5i=1 f(xi)
• . (b) f(x1, . . . , x5) = (c1,2c2,3c3,4c4,5)(c1,3|2c2,4|3c3,5|4)(c1,4|2,3c2,5|3,4)(c4,5|1,2,3) ∏5i=1 f(xi)

(c) f(x1, . . . , x5) = (c1,2c1,3c3,4c3,5)(c2,3|1c1,4|3c4,5|3)(c2,4|1,3c1,5|3,4)(c2,5|1,3,4) ∏5i=1 f(xi) Figure 2: Three examples of 5-dimensional copula constructions: (a) a D-vine; (b) a C-vine with the condition on T1 that u1 connects with every uj for j ̸= 1; The n-copula density is written below each graphical model.

and second using copulas due to Sklar's theorem (7):

From these two, we can derive that

by canceling out the first term f1(x1) in both. We can generalize this as

For the remaining part of this section variables will be omitted in probability densities and distributions since their subscripts give all the information. Thus, the last expression can be rewritten as

Next, let us consider three variables x1, x2 and x3. Like before, there are several factorizations of the joint density, f123, due to the chain rule. For instance, it could be factorized as: f1f23|1 or

Expanding f23|1 using Sklar's theorem we get

where c23|1 denotes the copula density c(F2(x2|x1), F3(x3|x1)).

According to (11) f3|1 = f3c31, hence by replacing f(3|1) in (13) we get:

Since (14) and (12) are both factorizations of the joint f123, by equating (12) and right hand side of (14) and canceling terms on both sides we get:

It is convenient to generalize this last expression, (15), as:

Finally consider a last example with four variables and the two following factorizations:

.

Sklar's theorem applied to f34|12, using (16), and straightforward manipulation leads to

Again, the last expression can be generalised as:

Hence, using (11), (16) and (17), any factorization of the 4 variable joint density can be expressed in terms of bivariate copulas. For instance:

A vine is a graphical representation of one factorization of the n-variate probability distribution in terms of n(n - 1)/2 bivariate copulas by means of the chain rule. It consists of a sequence of levels and as many levels as variables. Each level consists of a tree (no isolated nodes and no loops) satisfying that if it has n nodes there must be n-1 edges. Each node in tree T1 (level 1 ) is a variable and edges are couplings of variables constructed with bivariate copulas. Each node in tree T2 (level 2) is a coupling in T1, expressed by the copula of the variables; while edges are couplings between two vertices that must have one variable in common, becoming a conditioning variable in the bivariate copula. Thus, every level has one node less than the former. Once all the trees are drawn, the factorization is the product of all the nodes.

An example with 5 variables is given in Figure 2; showing two different vines together with their resulting joint pdf factorization.

For the sake of clarity we give the details of how to construct the panel of the left and how to derive its expression.

1. Construct T1, the first tree, in two steps:
2. (a) Draw in a row one node for each variable.
3. (b) Draw the edges by linking two adjacent nodes with a copula.

2. Construct T2 as follows:

• (a) Draw in a row one node for each edge in T1. Write the variables of the edges in T1 inside the nodes. For example in Figure 2 (left), the first edge in T1 is c1,2 which makes the first node of T2 1, 2.

• (b) Draw the edges by linking two adjacent nodes with a copula.

Write the bivariate copula such that the variables are those that only appear in one node or another, conditioning to the variables that are repeated in both; e.g. if nodes are C(1, 2) and C(2, 3), the conditioning variable is x2 because only 2 appears in both.

1. Repeat until the last tree, which has only one node.
2. The factorization is the product of every node.

□

Vines were initially presented in [Bedford and Cooke, 2001], [Cooke et al., 2007]. A comprehensive compilation of vine methods can be found in [Kurowicka and Joe, 2011]. In a nutshell, each vertex in Tj with j > 1 is a coupling in Tj-1 expressed by a bivariate copula with j - 2 conditioning variables which are those in common in the vertices coupled at Tj-1. Edges in Tj are formed between 2 vertices of Ti that add another common variable. The factorization is just the product of all the edges in all the trees, times the product of all the marginals. Let cj(v1, v2|w) be the density of an edge in Tj, so that v1, v2 ∈ {u1, . . . , ud} and v1 ̸= v2 are the variables of that edge, conditioned by the set w = {w1, . . . , wj-1} ∈ {u1, . . . , ud}, with wk ̸= v1 ̸= v2 for 1 ≤ k ≤ j - 1. Then, if E is the set of all edges in the vine, the factorization of the graphical model is:

In practice, however, it is usually recommended to avoid constructing all the trees because a full vine will only represent the actual underlying joint pdf if the bivariate copulas used are rightly chosen and accurately estimated [SalinasGuti´errez et al., 2010],[Haff et al., 2010].

In this paper we use two vine constructions: Drawablevines (D-vines) and Canonical-vines (C-vines), which are constructed according to different rules. All the models were fully obtained, in the sense that all the trees are developed, so that a study of how the depth influences the prediction can be done. The methodology in all of them is similar:

1. Transform all the variables by means of their marginals. In other words, compute

and compose the matrix u = u1, . . . , un, where ui are their columns.

1. For the construction of the first tree T1, assign one node to each variable and then couple them by maximizing the measure of association considered. Different vines impose different constraints on this construction. When those are applied different trees are achieved at this level. Figure 2 shows that D-Vine and C-Vine lead to different T1.
2. Select the copula that best fits to the pair of variables coupled by each edge in T1. Details about how this step is carried out are given in next subsection.
3. Let Cij(ui, uj) be the copula for a given edge (ui, uj) in T1. Then for every edge in T1, compute either v1 j|i = ∂/∂ujCij(ui, uj) or similarly v1 i|j , which are conditional cdfs. When finished with all the edges, construct the new matrix with v1 that has one less column u.
4. Set k = 2.
5. Assign one node of Tk to each edge of Tk-1. The structure of Tk-1 imposes a set of constraints on which edges of Tk are realizable. Hence the next step is to get a linked list of the accesible nodes for every node in Tk.
6. As in step 2, nodes of Tk are coupled maximizing the measure of association considered and satisfying the constraints impose by the kind of vine employed plus the set of constraints imposed by tree Tk-1.
7. Select the copula that best fit to each edge created in Tk.
8. Recompute matrix vk as in step 4, but taking Tk and vk-1 instead of T1 and u.
9. Set k = k + 1 and repeat from 6 until all the trees are constructed.

The rest of this section is devoted to the selection of the copula and the particular constraints for each vine.

Copula function selection

In step 4 above, there are multiple options for choosing the copula between a pair of variables. In this paper we consider three parametric copulas 1 : Clayton, Frank and Gumbel because with these we cover a wide range of tail dependences.There is a number of ways to evaluate which one of the copulas fits better to a dataset. Two of the most popular methods are to compare the empirical density function of the copula with the theoretical one [Genest and Rivest, 1993], and to compare the upper or lower tail functions [Venter, 2001].

In this paper we employ a triple-check fitting based on the latter. For this, we first compute an empirical Copula (nonparametric). We then numerically compute upper and lower tails given this empirical Copula and calculate the area under the tails given by al and au. The upper and lower tail concentration functions are respectively defined as:

Hence, given out three candidate copulas (Clayton, Frank and Gumbel), the procedure proposed for selecting the one that best fit to a dataset of pairs {(uj , vj )}j=1,2,..., is as follows:

1. Estimate the most likely parameter θ of each copula candidate for the given dataset.
2. Construct R(z|θ). Calculate the area under the tail for each of the copula candidates.

1 These are three popular archimedean families that can be found in many mathematical packages such as R or Matlab.

1. Compare the areas: au achieved using empirical copula against the ones achieved for the copula candidates. Score the outcome of the comparison from 3 downto 1, 3 being the best and 1 is the worst.
2. Proceed as in steps 2- 3 with the lower tail and function L.
3. Finally the sum of empirical upper and lower tail functions is compared against R + L. Scores of the three comparisons are summed and the candidate with the highest value is selected.

D-Vine Model.

In a D-Vine every node in tree T1 has degree 2 except two nodes, with degree 1, which can be seen as the extremes. Figure 2a shows an example of D-Vine of five variables. The advantage of D-Vine is that, once T1 is constructed, it uniquely determines the rest of the trees that compose the vine. Hence, learning the model is the task of finding the best assignment of variables to nodes in T1. Throughout this paper, Kendall's τ is employed as measure of association to decide how to couple nodes. The procedure is then as follows.

Let u be the matrix of the marginal cdf values of the dataset x;

1. Compute Mτ = [τij], with i > j, the matrix of Kendall's τ between every possible coupling of variables. Notice that Mτ only requires values in its upper triangular part because it is symmetric.
2. Find ˆ τ = max ([τij]). Let (left, right) be the coorτ .
3. dinates of ˆ Then initialize T1 = [left, right].
4. For k = 1 to m-1
5. (a) Set left = T1(1) and right = T1(end)
6. (b) Set the leftth and rightth columns of Mτ to zero.
7. (c) Find leftNew, the variable that couples with left with maximum Kendall's τ = τL. Similarly, find rightNew, the variable that couples with right with maximum Kendall's τ = τR.
8. (d) If τL > τR, then T1 = [leftNew, T1]; else T1 = [T1, rightNew]

C-Vine Model

In a C-Vine, for every tree, one anchor node is connected with all the others. Figure 2b shows an example with five variables. If the anchor node of T1 is the variable of the site of interest, then edges represent its dependence with respect to the rest of variables. The criterion for selecting the anchor node of Tk, for k > 1 is the following:

1. Set i = 1 and compute τij, the Kendall's of the data associated to node i and node j in Tk, for every j ̸= i.
2. Do τ[i] = ∑j τij
3. Repeat for all i.
4. The anchor node = arg max (τ[i]).

i

5 Results and discussion

In this section, we present the results obtained using the described wind resource assessment techniques on data acquired from the roof top anemometers at the Boston Museum of Science. We also examine the improvement in performance of each of the algorithms as more data is made available in the form of 3, 6 and 8 months of training sets. Additionally we study how much benefit is obtained when more expressive copula models are employed. To this end, we constructed the following 5 multivariate copulas and multiple regression (LRR).

Multivariate copulas constructed

Results are presented in Tables 1a-c for D3, D6, D8 respectively. Each table shows the KL distance between the ground truth distribution and the distribution estimated based on the predictions provided by each technique for the year 2 dataset (the test dataset) and for every bin. Their right-most column is the sum of the KL distance per bin, which gives an overall performance measure for every model. Tables are sorted according to this value. In addition, the minimum value of each bin and the model that attains most of these minimums have been highlighted.

5.1 Comparison of algorithms

First we compare algorithms when the same amount of data is available to each one of them for modeling. It is clear from Tables 1a-c that model LRR is the worst one with a large margin between it and the second and third worst, which are nGaussian and 4-t D-Vine. The remaining models attain very good results, with KL distances that range from 0.02 to 0.4. This result suggests that the model needs to incorporate a variety of dependence structures. Figure 4's heatmaps make it easier to extract qualitative conclusions across bins, models and training set size. Bands corresponding to C-vine models are, globally, darker than the rest for D3, D6 and D8. The figure indicates that models are better for west oriented bins (7-12) than for east oriented ones (1-6).

For a combined comparison, we consider two metrics for every model: (i) the sum of KL distance for every bin (rightmost column in Tables 1a-c), and (ii) the number of bins with minimum KL distance vs other models (highlighted cells). Sometimes the incomplete version of a vine performs better, or at least, equal to the full version. The most plausible explanation for this is that every tree added to a construction requires new copula estimations which always introduces another source of possible errors. In the C-Vine model the variable at the site of interest is the anchor node of the first tree. It influences every node in the whole C-vine which seems to explain the model's good performance compared with the incomplete D-Vine model. Because this variable is much less influential in D-Vine, one can expect to need more depth in order to attain similar results, as indeed happens.

As additional result, Figure 3 shows an example of truncated C-Vine; the one that corresponds to 8 month training

KL distance with 3 month training set

(a)

KL distance with 6 month training set

(b)

KL distance with 8 month training set

(c)

Figure 3: C-Vine truncated in the 4th level constructed out of the data from the 8 month training set and the 1st directional bin. Numbers represent the fourteen airports, and Y is the site of interest (Museum of Science, Boston). The bold edge in every tree represents the couple with highest Kendall's τ , and τ decreases clockwise.

set for the first directional bin; which is highlighted in Table 1c as its best model. The bond between nodes with highest Kendall's tau is represented with a bold edge and edges are deployed clockwise as τ decreases.

5.2 Increasing the data available for modeling

We next examine the robustness of each technique as progressively more data is made available to it for modeling. Figure 5a compares the sum of KL distances of each model for training sets D3, D6 and D8. C-Vine does not significantly change as more/less data is incorporated whereas the other models get worse with more data. This may indicate overfitting which could be a disadvantage of such sensitive, tunable models.

When we examine the minimum KL distance attained with each train set for every bin (not shown),the highest difference in KL distance is less than 0.1. In other words, there is always at least one model out that performs similar to the best one in case of lost data or when more data is available.

6 Conclusions

In this paper we presented copula based approaches for Wind resource estimation. Copula based approach allow us form a joint distribution with Weibull marginals and allow us capture non-linear correlations between the variables. In addition, we presented a methodology to construct a variety of copula models by factorizing the joint in different ways. With its ability to capture long tails and tail dependencies these models allowed us to estimate the wind resource at the new site with as little as 3 months of data. This is a significant achievement in the wind resource estimation domain where ability to estimate the wind resource accurately in less amount of time allows better planning. Such estimation from reduced amount of time/data is highly beneficial for offshore wind technology development where site-measurement campaigns are extremely expensive.