2.1. A model based on Gaussian process regression

To model the normal spatiotemporal pattern across a wind farm, a flexible regression algorithm is required. In this work, the model of choice is Gaussian process regression (GPR).

Over the last decade, GPR-based methods have gained huge popularity in the SHM community (Cross, Reference Cross2012; Papatheou et al., Reference Papatheou, Dervilis, Maguire, Antoniadou and Worden2015; Bull et al., Reference Bull, Gardner, Gosliga, Maguire, Campos, Rogers, Haywood-Alexander, Dervilis, Cross, Worden and Mao2020; Lin et al., Reference Lin, Worden, Maguire, Dilworth and Mains2020b; Rogers et al., Reference Rogers, Gardner, Dervilis, Worden, Maguire, Papatheou and Cross2020). As a probabilistic model, a GPR is able to provide a predictive distribution, which indicates both the mean predicted values and their confidence interval (Rasmussen and Williams, Reference Rasmussen and Williams2006). In the context of PBSHM, the predictive distribution also accounts for the population variance, that is, the difference between individual structures within a population that arises from benign sources, such as turbulent boundary conditions and manufacturing tolerances, so that a general form can be established across multiple structures of the same “type” (Bull et al., Reference Bull, Gardner, Gosliga, Maguire, Campos, Rogers, Haywood-Alexander, Dervilis, Cross, Worden and Mao2020). Another advantage of GPR lies in the fact that it is nonparametric and well suited for learning complex relationships, as opposed to a parametric model whose complexity may be restricted by a predefined functional form.

As an extension to multivariate Gaussian distributions that describe random variables of finite dimensions, a GPR is held over functions, with the function values $ f(X) $ specified at a potentially infinite number of inputs $ X $ . The fact that a GPR can provide consistent inferences over any finite number of points makes it a powerful tool to learn functional relationships (Rasmussen and Williams, Reference Rasmussen and Williams2006). A GPR is specified by a mean and a covariance function:

(1) $$ f(X) \sim \mathcal{GP}\left(m(X),k\Big(X,{X}^{\prime}\Big)\right). $$

Conventionally, the mean function is assumed to be zero to simplify the formulation, and the covariance function takes the form of squared exponential (SE), which is the most commonly used kernel in describing smooth data variations:

(2) $$ {k}_{\mathrm{SE}}\left(X,{X}^{\prime}\right)\hskip0.35em =\hskip0.35em {\sigma}_f^2\;\exp \left(-\frac{1}{2{l}^2}{\left|X-{X}^{\prime}\right|}^2\right). $$

The process of model training involves the optimization of hyperparameters, which, according to Eq. (2), include process variance $ {\sigma}_f^2 $ and input length scale $ l $ . Note that characteristic length scales are used in the current model, meaning that each input dimension is assigned a different length scale; for input data with a dimension $ D $ , there is a total of $ D $ length scales $ {l}_1,\dots, {l}_D $ . These hyperparameters are optimized here by maximizing the marginal likelihood via gradient-based optimization methods (Rasmussen and Williams, Reference Rasmussen and Williams2006).

With the optimized hyperparameters, a predictive distribution corresponding to the unseen inputs, $ {X}^{\ast } $ , can be given by specifying the predictive mean vector, $ {\overline{\mathbf{f}}}^{\ast } $ , and covariance matrix, $ \operatorname{cov}\left({\mathbf{f}}^{\ast}\right) $ :

(3) $$ {\overline{\mathbf{f}}}^{\ast}\hskip0.35em =\hskip0.35em K\left({X}^{\ast },X\right){\left[K\left(X,X\right)+{\sigma}_n^2I\right]}^{-1}\mathbf{y}, $$

(4) $$ \operatorname{cov}\left({\mathbf{f}}^{\ast}\right)\hskip0.35em =\hskip0.35em K\left({X}^{\ast },{X}^{\ast}\right)-K\left({X}^{\ast },X\right){\left[K\left(X,X\right)+{\sigma}_n^2I\right]}^{-1}K\left(X,{X}^{\ast}\right), $$

where $ K\left(\cdot, \cdot \right) $ denotes a covariance matrix obtained from Eq. (2), and $ {\sigma}_n^2 $ is the noise variance associated with the noisy observations $ \mathbf{y} $ . Since the goal here is to predict the noisy target $ {\mathbf{y}}^{\ast } $ , the covariance for noisy predictions should be $ \operatorname{cov}\left({\hat{\mathbf{y}}}^{\ast}\right)\hskip0.35em =\hskip0.35em \operatorname{cov}\left({\mathbf{f}}^{\ast}\right)+{\sigma}_n^2I. $

For more details on GPR formulation/implementation, readers are referred to Rasmussen and Williams (Reference Rasmussen and Williams2006).

2.2. Functional relationships captured by the model

As mentioned earlier, the mapping method aims to capture the spatiotemporal correlation between the EOVs and turbine performance. The choices of features are given as follows: power is chosen as the target feature, since it is a direct result of environmental variations as well as operational control, and, consequently, a direct measurement of turbine performance; features representing the EOVs include the spatial coordinates at all turbine locations and the wind speed at a subset of locations. A GPR model is trained to predict the power time series at all turbine locations in a wind farm, from the corresponding spatial coordinates and wind speed at a subset of locations. In other words, the model intends to capture the wind–power relationship and how it is affected by the wake-induced correlation across space and time.

The wind–power correlation can be summarized by a power curve, with an example shown in Figure 1. It is seen that power is strongly correlated with wind speed, forming an approximately sigmoidal (or cubic within a specific range) relationship. The shape of the curve is piecewise as a result of the turbine control strategy; thus, the modeling of a full power curve is a complicated task in itself. The use of a simple GPR model, without incorporating physical knowledge or partitioning the input space, might result in reduced accuracy in regions of control intervention (e.g., above-rated power; Rogers et al., Reference Rogers, Gardner, Dervilis, Worden, Maguire, Papatheou and Cross2020). Therefore, it is unlikely that a simple GPR is able to capture the full wind–power correlation as well as the spatiotemporal correlation induced by wakes. As a result, the data used in the following analysis belong to a restricted set of operating conditions, in which the wind speed variations are kept within the range in-between the cut-in and rated values (Figure 1). An extension to modeling the full wind–power correlation may be useful once the current model is proved valid.

Figure 1. An example of a power curve.

At the risk of repetition, it should be clarified that power-curve monitoring is not part of the aim of the mapping method (and the turbine-by-turbine wind–power correlation is not modeled directly). Instead, wind speed and power are the chosen features that represent EOVs and structural response, respectively. The aim of the mapping method is to achieve EOV modeling in a spatiotemporal setting, and the aim of this paper is to demonstrate the potential of this method as a performance indicator for a wind farm.

The GPR model not only captures how power changes with wind speed temporally, but how this correlation varies spatially. The spatial aspect comes from the wake effect in the wind farm, which affects wind speed and power in a similar manner. As illustrated in Figure 2, the maximum power output in a farm is found at the front row of turbines; the power levels decrease as the wind progresses further downstream across the farm. This spatial correlation is modeled by using as input variables (a) the spatial coordinates of all turbine locations and (b) the wind speed time series at a fixed set of reference locations. The reference locations are selected to ensure that the model learns the spatial and temporal correlations through interpolation (i.e., to avoid extrapolation).

Figure 2. Example power variations across the Lillgrund wind farm corresponding to different wind directions, which are indicated by the blue arrows.

Figure 2 also shows that the spatial pattern is sensitive to wind direction, which changes with time. This time-varying spatial pattern is also captured by the model.

To conclude, the GPR model captures (a) the correlation between wind speed at reference locations and power at each output turbine, (b) how the wind–power correlation varies spatially, and (c) how the spatial pattern changes with wind direction (which is time-varying).

2.3. Discrepancy measures and thresholds

After establishing a model of normal condition, the expected normal turbine performance during a test period can be predicted, and the discrepancy between the prediction and the measurement is computed. As a population-based method, the mapping method detects anomalies via comparisons across the structures in a population. Therefore, for each turbine, the difference between the predicted and measured power time series needs to be summarized into a single discrepancy value to facilitate the comparison between turbines. Two discrepancy measures are used for the purpose of this study, namely normalized mean-squared error (NMSE) and mean standardized log loss (MSLL), which are error metrics commonly used for model quality evaluation in the literature of GPR.

The NMSE is given by

(5) $$ {\mathrm{NMSE}}_s\hskip0.35em =\hskip0.35em \frac{100}{T{\sigma}_{y^{\ast}}^2}\sum \limits_{t\hskip0.35em =\hskip0.35em 1}^T{\left({y}_{t,s}^{\ast }-{\overline{\mathrm{f}}}_{t,s}^{\ast}\right)}^2, $$

where $ y $ and $ \overline{\mathrm{f}} $ denote the measured and predicted (mean) target variables, respectively. The superscript $ {}^{\ast } $ indicates the test dataset, and the subscripts $ t $ and $ s $ denote the time instance $ t\hskip0.35em =\hskip0.35em 1,\dots, T $ and the turbine index $ s\hskip0.35em =\hskip0.35em 1,\dots, S $ . The mean-squared error term is normalized by the term $ 100/{\sigma}_{y^{\ast}}^2 $ . $ {\sigma}_{y^{\ast}}^2 $ represents the variance of the measured target variable in the test set across all turbines, in order to remove error sensitivity to the scale of testing data. The scaling factor 100 converts the error values into percentages. A score of 100% means that the predictions are equivalent to the target mean of the test data. The threshold of NMSE is therefore set to be 100%, below which the model provides better predictions than the mean of test data and is assumed to have captured some correlations between input and output. Note that the calculation of NMSE does not include the information about GPR predictive variance.

The second metric, MSLL, includes both the mean and variance of the GPR predictive distribution in its formulation (Rasmussen and Williams, Reference Rasmussen and Williams2006):

(6) $$ {\mathrm{MSLL}}_s\hskip0.35em =\hskip0.35em \frac{1}{2T}\sum \limits_{t\hskip0.35em =\hskip0.35em 1}^T\left(\log \left(2{\pi \sigma}_{{\hat{y}}_{t,s}^{\ast}}^2\right)+\frac{{\left({y}_{t,s}^{\ast }-{\overline{\mathrm{f}}}_{t,s}^{\ast}\right)}^2}{\sigma_{{\hat{y}}_{t,s}^{\ast}}^2}\right)-\frac{1}{2T}\sum \limits_{t\hskip0.35em =\hskip0.35em 1}^T\left(\log \left(2{\pi \sigma}_y^2\right)+\frac{{\left({y}_{t,s}^{\ast }-\overline{y}\right)}^2}{\sigma_y^2}\right). $$

The first term in Eq. (6) is the negative log probability that indicates how likely the target measurements are to be predicted using the trained model. Here, the GPR-predicted variance for a noisy target is $ {\sigma}_{{\hat{y}}_{t,s}^{\ast}}^2\hskip0.35em =\hskip0.35em \operatorname{diag}\left(\operatorname{cov}\left({\hat{\mathbf{y}}}^{\ast}\right)\right)\hskip0.35em =\hskip0.35em \operatorname{diag}\left(\operatorname{cov}\left({\mathbf{f}}^{\ast}\right)+{\sigma}_n^2I\right) $ . The second term in Eq. (6) represents the (negative log) probability obtained by treating the training mean $ \overline{y} $ and variance $ {\sigma}_y^2 $ as the model. It standardizes the first term such that $ \mathrm{MSLL}>0 $ if the GPR predictions are worse than the trivial model of training mean and variance. The threshold of MSLL is set to be 0, below which the model prediction is considered acceptable.

In the analysis that follows, it is demonstrated that the abovementioned discrepancy measures, NMSE and MSLL, and their corresponding thresholds can be used for anomaly detection.

3.1. Datasets

The SCADA data are collected from an offshore wind farm, Lillgrund. In this study, time series data refer to a series of mean values summarizing every 10-minute section of the raw data streams. Detailed description and analysis of the wind farm can be found in Jeppsson et al. (Reference Jeppsson, Larsen and Larsson2008) and Dahlberg (Reference Dahlberg2009).

A GPR-based model of normal condition is trained to predict power from spatial coordinates and wind speed, as mentioned in Section 2.2. One of the inputs to the model is the spatial locations, which are given as the Cartesian coordinates of each turbine position. Another input feature is the wind speed at a number of reference turbines. A set of 10 reference turbines is randomly selected across the wind farm, by manually choosing the minimum distance among them such that they are (roughly) evenly distributed across the entire space. For the purpose of this preliminary study, the only criterion behind the selection is to (visually) ensure that the reference locations are spread across the wind farm in order to avoid extrapolation in space, since a standard/zero-mean GPR model tends to predict less accurately during extrapolation (Cross and Rogers, Reference Cross and Rogers2021). More advanced methods for reference selection, such as Latin hypercube sampling, may be considered for future studies. The chosen reference locations are indicated by the rectangles in Figure 3.

Figure 3. Maps of (a) NMSE and (b) MSLL averaged across a testing data period (of 2 hours 40 minutes) that describes normal operational conditions. The incoming wind direction is indicated by the blue arrow, and the reference turbine locations are marked with the blue rectangles. The turbine numbers (1–48) are also indicated, outside the circles.

The training and testing datasets are selected based on the following criteria. First, the data on wind speed and power used in the case study are confined to the range between cut-in and rated values (Figure 1), for the reasons mentioned in Section 2.2. Second, the GPR model includes the effect of multiple wind directions as part of the definition of normal condition; thus, the training data cover a range of wind directions. It is noted that the training data are selected so as to avoid rapid/extreme fluctuations in wind directions. If wind direction changes too rapidly, the turbines may not be given enough time to respond to the change, and the resulting spatial pattern may not be reflective of the wind direction. Given that the mapping method aims to capture how spatial patterns change with wind directions, data describing steady changes in wind directions are necessary for model training. A summary of the training and testing datasets used in the case study is given in Table 1. Note that the numbers of data points for the training period in Table 1 correspond to the numbers obtained after subsampling the original SCADA data by a factor of 2 (in order to reduce computation time), whereas the testing data are not subsampled. Thus, to make test predictions, the model interpolates not only in space but in time. The training data correspond to about 25.5 hours’ worth of data for all 48 turbines in the Lillgrund wind farm, and the test datasets 1 and 2 correspond to 2.8 and 2.2 hours’ worth of data, respectively. Note that the individual data sections are the longest continuous periods during which the mentioned criteria (of wind speed, power, and wind direction) are satisfied.

Table 1. A summary of training and testing datasets.

3.2. Capability to predict normal power production

The trained GPR model is used to predict the power output during a testing period, when all turbines are operating under normal conditions. This example refers to Test Set 1 in Table 1. The differences between measurements and predictions are summarized by the maps of NMSE and MSLL shown in Figure 3. It is seen that no specific spatial pattern can be found in the NMSE map (Figure 3a), whereas the MSLL results tend to get better toward the downstream side of the farm (Figure 3b). Such a spatial trend in the MSLL values is arguably caused by the calculation of the MSLL requiring an indirect comparison of the mean and variance between training data and GPR predictions; that is, the two terms in Eq. (6) correspond to the log loss obtained by comparing the testing data with GPR predictions and training statistics, respectively. Since the predicted mean power levels tend to decrease in the direction of the wind, the MSLL values change in the wind direction as well, but the way the MSLL changes depends on how far away the predictions are from the training mean. Note that the confidence levels of model predictions remain roughly constant throughout the farm. Overall, the error values are well below the corresponding thresholds (100% for NMSE; 0 for MSLL). At reference turbine locations, that is, where wind speed measurements are provided as the model inputs, there is a tendency to obtain relatively better NMSE results, while whether reference inputs are provided does not seem to affect the MSLL values. The difference between the two error metrics mainly stems from the fact that only the mean predictions are used in the NMSE, whereas both the predictive mean and variance are included in the calculation of the MSLL.

Examples of time series predictions are visualized in Figure 4, where the GPR prediction constitutes a mean and a confidence interval (given as three times the GPR predictive standard deviation $ {\sigma}_{{\hat{y}}^{\ast }} $ ). Figure 4a illustrates one of the best predictions, excluding those at reference turbines. It is seen that the mean prediction follows the measured trend almost exactly, even at regions of relatively wider confidence intervals (from Data Point 11 onward). Figure 4b demonstrates the model prediction with a medium level of errors, where there are periods of mismatch between the mean prediction and the measurement, but the amount of deviation is considered small compared with the confidence interval. The largest-error prediction is illustrated in Figure 4c; it can be seen that the mean prediction deviates from the measurement throughout the time window while still loosely following the trend. In regions of large deviation (Data Points 5–9), the measured data are close to the boundary of the confidence interval, which leads to the relatively large MSLL as well as NMSE. In summary, the model is able to provide mean predictions that follow (at least roughly) the measured trends, as well as confidence intervals that capture all the benign power variations in the given testing dataset.

Figure 4. Examples of the predicted and measured power time histories for the normal testing set, in cases of (a) small, (b) medium, and (c) large errors.

3.3. Performance indicator for wind farms

Having established an appropriate model of normality, the model is used to predict the power production in a test set that contains potential performance anomalies. Figure 5 shows the discrepancy between predictions and measurements in terms of maps of NMSE and MSLL. The color schemes of the error maps are designed such that the turbines with an error value above the threshold are highlighted (in warm colors). The thresholds for NMSE and MSLL are 100% and 0, respectively, with reasons given in Section 2.3. It is seen that five turbines are highlighted as candidate anomalies in the NMSE map, that is, Turbines 22, 29, 33, 38, and 48, and there is one additional candidate anomaly, Turbine 25, highlighted in the MSLL map. Among these highlighted turbines, Turbine 48 is disqualified from being a potential anomaly, since it was manually switched off during the time period being considered. This powers down results in extremely high values of NMSE and MSLL. Although this is a known control action, it provides very basic validation that the model is able to detect anomalous behavior. The other highlighted turbines will now be investigated. Note that the dataset available is currently unlabeled; that is, the true reason for suspected anomalies is unknown—here some possible reasons for the observed behavior will be explored.

Figure 5. Maps of (a) NMSE and (b) MSLL averaged across a 2-hour testing period with potential anomalies. The incoming wind direction is indicated by the blue arrow, and the reference turbine locations are marked with the blue rectangles. The turbine numbers (1–48) are also indicated, outside the circles.

The time histories of the predicted and measured power for the candidate anomalies (Turbines 22, 25, 29, 33, and 38) are shown in Figure 6. For all candidates, there are obvious discrepancies between mean predictions and measurements throughout the time window. However, the trends of measured power are still captured by the confidence intervals most of the time. Interestingly, in all five cases, the model underpredicts the power produced, implying that there may be a factor that causes many turbines to produce more power than usual during this period of time. With respect to the error values, the further away the mean predictions are from the measurements, the higher the NMSE values, whereas higher MSLL values tend to result from scenarios when the measured data trends are consistently close to or beyond the uncertainty boundaries. In general, none of the candidate anomalies exhibits large fluctuations in power, and the relatively high errors seem to result from predictions with shifted means.

Figure 6. Time histories of the predicted and measured power for five highlighted candidate anomalies.

One reason why the actual power production deviates from model prediction is that the testing dataset describes a spatial correlation different from that depicted in the training data. This difference in spatial correlation is the main cause of the candidate anomalies at Turbines 25 and 38. In Figure 7, the training and testing wind–power correlations for the first three turbine rows, that is, the neighborhood of Turbines 25 and 38, are illustrated, with comparisons made available between training and testing data. A comparison between Figure 7a,b illustrates how the spatial correlations between Turbine 25 and its neighbors change from training to testing data. It is seen that, although Turbine 25 tends to generate relatively high power, given the same wind speed compared to its neighbors during training (Figure 7a), the amount of extra power it generates in the testing period becomes more distinct from the neighborhood (Figure 7b). However, what is shown in Figure 7b can (arguably) be merely a demonstration of the population variance; perhaps that is why Turbine 25 is flagged as a candidate anomaly by only one of the error metrics. For Turbine 38, the difference between training and testing data lies mainly in its relationship with the upstream Turbine 37. During training, the effect of wake shielding is apparent as the wind and power at Turbine 38 are significantly less than those at Turbine 37 (Figure 7c), which is not seen in the testing period (Figure 7d). Wake deviations due to yaw misalignment may be a reason for such difference, which remains hypothetical since the relevant data are unavailable to the authors. The difference in spatial pattern around Turbine 38 is captured by both error metrics.

Figure 7. Wind-power correlations associated with (a)-(b) Turbine 25 and (c)-(d) Turbine 38. The “Neighbors” are referred to turbines in the neighborhood of Turbines 25 and 38, in this case the first three rows of turbines in the wind direction (excluding other candidate anomalies).

Another cause of deviation from normal condition is an unexpected local nacelle direction. In particular, if a turbine turns away from the wake(s) of the upstream turbine(s), the power produced by the turbine may be underestimated by the model. Such an unexpected nacelle direction is found at Turbine 33. Normally, Turbine 33 should be affected by the combined wakes of the two upstream turbines (Turbines 31 and 32). However, during the testing period, this turbine is constantly facing Turbine 37, a region likely to be affected by fewer wakes. Therefore, Turbine 33 ends up with more power output than predicted, as seen in Figure 6d, and this overprediction has been successfully highlighted by both maps in MSLL and NMSE.

The remaining candidate anomalies, Turbines 22 and 29, stand out from their neighbors, as they experience unexpectedly high wind speed and, thus, produce more power. Figure 8 shows the wind–power correlations of Turbines 22 and 29 in relation to those of the neighboring turbines. As part of the normal spatial pattern acquired by the model, Turbines 22 and 29 generate relatively low power in comparison to their surrounding turbines (Figure 8a). In the testing period, the wind and power at the two turbines extend far beyond the ranges covered by their neighbors while maintaining the gradient of the wind–power correlations (Figure 8b). Therefore, the difference between training and testing data is likely to stem from changes in the environment rather than changes in the turbine systems. Upon further examination of the data under similar environmental and operational conditions, the unexpected high speeds at Turbines 22 and 29 are unlikely to be caused by unexpected nacelle direction or ambient temperature, nor are there seasonal or diurnal patterns in the occurrence of such phenomena. Given the existence of a central gap in the wind farm of interest (see Figure 5), more wind eddies may be generated around the gap region, bringing about more short-term turbulence that causes the occasional high wind around Turbines 22 and 29. The fact that these infrequent environmental variations are not described by the training data has led to the seemingly sensitive detection results. However, the occasional exposure to wind eddies that give rise to counterintuitive wake patterns may subject Turbines 22 and 29 to higher fatigue loads and, thus, a higher risk of anomalous performance. The model results, therefore, provide insight into this potential risk.

Figure 8. Wind-power correlation for a subset of turbines toward the back of the farm, that is, Turbines 13–15, 21–23, 28–30, 35, and 36, obtained from (a) training and (b) testing data.

3.4. Effect of reference wind speed inputs

As previously stated, the mapping model is designed to predict the power across a wind farm, as a function of turbine locations and wind speed at a fixed set of reference turbines. In the previous case study, 10 reference turbines are chosen, with their locations indicated in Figures 3 and 5. This section focuses on investigating how the numbers of reference turbines may affect model accuracy. For a given number of reference turbines, the locations were randomly selected five times, and a GPR was trained each time. The NMSE results corresponding to different numbers of reference turbines are illustrated in Figure 9. It is seen that, as the number of reference turbines increases, there is an approximately exponential decrease in both the NMSE values and the error variances. If the reference locations are to be optimized, the NMSE improvements will be larger at fewer numbers of references. Optimized reference locations are considered to bring significant NMSE reduction when the number of references is smaller than or equal to 8. The average NMSE level for 10 references is around 30%. For the set of reference turbines used in the case study, the overall NMSE is around 20% across all normal testing datasets. Hence, in the current study, an optimized set of reference locations may not bring significant improvements in predictive accuracy. Nonetheless, an optimization algorithm, such as a genetic algorithm (Worden et al., Reference Worden, Manson, Hilson and Pierce2008), will be necessary if fewer references are to be selected in future studies.

Figure 9. The NMSE for various numbers of reference turbines.