Gaussian process power curve models incorporating wind turbine operational variables

The IEC standard 61400 − 12 − 1 recommends a reliable and repeatable methodology called ‘binning’ for accurate computation of wind turbine power curves that recognise only the mean wind speed at hub height and the air density as relevant input parameters. However, several literature studies have suggested that power production from a wind turbine also depends significantly on several operational variables (such as rotor speed and blade pitch angle) and incorporating these could improve overall accuracy and fault detection capabilities. In this study, a Gaussian Process (GP), a machine learning, data-driven approach, based power curve models that incorporates these operational variables are proposed in order to analyse these variables impact on GP models accuracy as well as uncertainty. This study is significant as it find out key variable that can improve GP based condition monitoring activities (e.g., early failure detection) without additional complexity and computational costs and thus, helps in maintenance decision making process. Historical 10-minute average supervisory control and data acquisition (SCADA) datasets obtained from variable pitch regulated wind turbines, are used to train and validate the proposed research effectiveness The results suggest that incorporating operational variables can improve the GP model accuracy and reduce uncertainty significantly in predicting a power curve. Furthermore, a comparative study shows that the impact of rotor speed on improving GP model accuracy is significant as compared to the blade pitch angle. Performance error metrics and uncertainty calculations are successfully applied to confirm all these conclusions.


Introduction
Over recent decades, wind power has experienced fast development and emerged as a viable and cost-effective alternative to conventional power generation. Reflecting the rapid installation of wind turbines, the European Wind Energy Association (EWEA), (Moccia et al., 2011), has projected that the total electricity demand in Europe could well reach 400 GW by 2030, of which 28.5% will be supplied by wind power. However, the overall cost of energy (CoE) from an offshore wind farm remains high with significant Operation and Maintenance (O&M) costs comprising a significant contribution to total costs and making offshore power financially less attractive than it would be otherwise. For example, recent research (Ioannou et al., 2018) revealed that O&M cost can make up 25% of the overall lifecycle cost of generation in the case of offshore wind turbines (WTs). The premature failure of critical components of WTs can cause underperformance and power production loss, as well as * Corresponding author.
E-mail address: ravi.pandit@strath.ac.uk (R.K. Pandit). significant downtime and reduction in availability (Scheu et al., 2019;Leimeister and Kolios, 2018;Scheu et al., 0000). Therefore, it is in the interest of wind farm owners and operators to detect failures in a timely manner and thus prevent catastrophic damage and so optimise maintenance and availability, making offshore wind a more profitable business.

Recent works on SCADA data based data-driven techniques
Supervisory control and data acquisition (SCADA) systems store data that give valuable information about the operational status of WTs at no additional cost. Therefore, SCADA-based condition monitoring is a cost effective approach and has attracted significant interest in recent years from researchers and industry practitioners, aiming to solve problems such as early failure detection and performance/condition monitoring, (Castellani et al., 2017;Martinez-Luengo et al., 2016). Researchers have used a range of statistical techniques for WT condition monitoring based on SCADA data; these can be classified as parametric and non-parametric methods. Parametric methods (e.g., polynomial regression) are usually restricted in their effectiveness by Nomenclature f (x) Gaussian Process distributed function for given data x K ( x, x ′ ) Covariance between x and x ′ data points m (x)

Calculated air density values V C
Corrected wind speed in m/sec V M Measured wind speed in m/sec their dependency on specific mathematical equations. In contrast, non-parametric methods refer to data-driven models that do not impose any pre-specified model and generally perform better than parametric models in terms of accuracy and ease of implementation. Kusiak and Zhang (2010) and Zhang and Kusiak (2012), used SCADA data of wind speed, blade pitch angle, and generator torque for detecting faults based on nonparametric models that include ANN (Artificial Neural Network), ANN ensemble, boosting regression trees, SVM, random forest, k-nearest-neighbour ANN, among others. Furthermore, Zhang and Wang (2014), proposed an ARX (Auto-Regressive with eXogenous input) ANN model to detect a main bearing fault, where nacelle temperature and turbine speed were used as exogenous inputs. In Lapira et al. (2012), three machine learning techniques (self-organising map, Gaussian mix model and ANN) are proposed and compared to evaluate the performance of wind turbines. Wang and Infield (2013), proposed a non-parametric technique called 'non-linear state estimation technique (NSET)' to model a healthy wind turbine gearbox using historical 10-minute average SCADA data. Welch's t-test was used in the fault detection algorithm, together with suitable time series pre-processing, for early detection of incipient anomalies in the turbine gearbox, i.e., before they reached a catastrophic stage. The constructed NSET fault detection algorithm was then compared with an ANN model and was found to have superior performance. Prior to Wang and Infield (2013), Guo (2012) proposed the NSET technique to model the generator bearing temperature but did not actually apply the technique to early failure detection. In Zhao et al. (2018), a deep auto-encoder (DAE) deep learning network based on SCADA data for component early anomaly detection and fault detection was proposed and its effectiveness verified by reported failure cases of wind turbine components (generator and gearbox). A brief but useful literature review of widely used machine learning/ data-driven approaches in WTs applications can be found in Stetco et al. (2018).

Works related to wind turbine power curve
The power curve is widely used to either monitor the performance of turbines or to provide wind power prediction for given wind speed, (De Giorgi et al., 2011). Therefore, accurate modelling of power curves is important for accurate forecasting of power from wind turbines and also for reducing O&M cost through early detection of changes in performance. Several techniques have been used for power curve modelling: parametric (Taslimi-Renani et al., 2016) and nonparametric (Manobel et al., 2018;Dongre and Pateriya, 2019), and also stochastic (Gottschall and Peinke, 2008) methods. Parametric models are defined by specific mathematical equations while nonparametric approaches do not enforce any prespecified conditions; as a consequence power curves based on nonparametric models closer to the measured power curve over a wide wind speed range. Popularly used nonparametric models that are extensively used for power curve modelling or related activities are GPs (Pandit and Infield, 2017), support vector machine , artificial neural networks (Marvuglia and Messineo, 2012) and Gaussian copula (Gill et al., 2012). Recent comprehensive studies on existing power curve monitoring techniques can be found in Marvuglia and Messineo (2012), Gill et al. (2012) and Astolfi et al. (2019). All the above research demonstrates the usefulness of machine learning and/or data-driven approaches to construct a feasible condition monitoring system for WTs out of which most of them used power curve knowledge to build a robust condition monitoring system, see for example, Castellani et al. (2017), Martinez-Luengo et al. (2016), Zhang and Wang (2014) and Lapira et al. (2012).

Works related to GP for WTs
A Gaussian Process (GP) is a data-driven, non-parametric, supervised machine learning technique used in regression as well as in classification problems. Recently, GPs are finding applications in solving wind turbines issues because of their ease of operation and flexibility of algorithm construction. For example, in Chen et al. (2014), authors describe a technique that includes a NWP (Numerical Weather Prediction) model, which is incorporated in a GPR (Gaussian Process Regression) model, to forecast wind speeds up to 1 day ahead. Then, corrected wind speeds are used to predict wind power using another GPR method. They used three SCADA datasets obtained from different wind farms located in China to validate the effectiveness of the proposed method. In Pandit and Infield (2018), GPR method based on SCADA proposed to detect a failure caused by yaw misalignment, which is then compared against the two binning based methods. The result shows that the GPR model was able to detect the failure with the alarm raised only 1.5 h after the fault occurred, while binning algorithms took 6 h and ∼4 h. Furthermore, in , the same authors proposed, and carried out performance comparisons for, three data-driven power curve models, namely, GP, SVM and Random Forest (RF) and found that the GP-based power curve model accuracy was better than for SVM and RF. In a different study, Li et al. (2019) proposed a SCADA data-based GPC (Gaussian Process Classification) model for WT diagnosis and then compared this with SVM. The results suggest that the GPC method is able to provide more accurate fault diagnosis results than the SVM on average. GP models have also been extended to solve power output control (Jin et al., 2014;Kou et al., 2015) and wind power forecasting (Yan et al., 2016;Lee and Baldick, 2014), for wind turbines. Despite promising results, GP applications related to the power curve have been limited to date and further work is needed to establish their potential for WTs condition monitoring and early failure detection related activities.

Scientific novelty and the importance of this research
In general, while constructing a power curve, all these datadriven techniques considered wind speed and air density as the relevant input parameters, following the IEC Standard 61400- 12-1 (2006), and ignored other relevant parameters such as blade pitch angle and rotor speed (called operational variables) that are known to affect power production. Although the latest version of the IEC standard, (IEC Standard 61400-12-1, 2017), highlighted the impacts of turbulence and wind shear and the author of Castellani et al. (2017) has analysed the impact of wind direction on power production, both have ignored the impact of the operational variables. Thus, it is unclear whether the inclusion of operational variables will improve power curve accuracy as well as reduce modelling uncertainty. Therefore, this paper explores a GP-based power curve model that incorporates these operational variables. Analysing the impact of operational variables on GP power curve and then finding most appropriate operational variable can not only improve the power curve but is also helpful for constructing robust GP power curve based fault detection algorithms. This is the motivation for the work presented here.
A framework for the proposed research is illustrated in Fig. 1 and described as follows. The SCADA data extracted from healthy turbines are first filtered and corrected for air density. The preprocessed data are then split into training and validation sets; the proposed GP model is fitted using the training datasets while validation datasets are used to validate the performance of the GP model. Here, operational variables (pitch angle & rotor speed) are incorporated along with air density corrected wind speed to train and validate the GP model. Performance comparison is the final stage of the proposed methodology. Statistical performance error metrics and uncertainty assessment are undertaken to answer the following research question: does the inclusion of operational variables improves the GP model accuracy, and if so, which among them is most significant in improving the accuracy?
The rest of this paper is organised as follows. Section 2 discusses WT power curve characteristics. Section 3 describes SCADA data description, pre-processing and air density correction. In Section 4, the GP methodology for power curve modelling is described. Incorporating operational variables and their impact on GP power curve model accuracy and uncertainty is investigated in Section 5, which has two subsections. Section 6 presents comparative performance studies for the GP model incorporating operational parameters. Finally, conclusion and future work related to this research is included in Section 7.

The characteristics of wind turbine power curve
The dynamic, nonlinear relationship between electrical power output and wind speed is defined by the power curve (see for example Fig. 2), and is in part governed by a cubic relationship between these variables. This relation can be expressed mathematically the following equation (Jin et al., 2014): , A is swept area (m 2 ) , C p is the power coefficient of WT and v is the hub wind speed (m/sec).
The power coefficient is a nonlinear function of the tip speed ratio (TSR) (λ) and blade pitch angle (β).
The electrical power output is not only dependent on wind speed but also influenced by other important variables such as atmospheric pressure, wind shear, rotor diameter, rotor height turbulence intensity, wind direction variability, both vertical and horizontal shear, atmospheric stability, drive train temperature and so on (IEC Standard 61400-12-1, 2006). The impact of these variables on performance can cause a mismatch between measured and manufactures' power curves. Though both show similar trends, with real data the measured power curve is more scattered. Furthermore, power curve-based techniques ignore vital information such as turbine subsystem temperature, environment status or turbine operational status. In addition to this, the power curve does not include technical details such as local terrain, wind direction, turbine wakes and other factors. However, in Castellani et al. (2017), IEC Standard 61400-12-1 (2017),  it is recognised that inclusion of these variables will significantly influence the power curve accuracy as well as uncertainty.
Power curve based fault detection techniques have been extensively applied and found to be effective, though most of the studies, such as Gill et al. (2012) and Pandit and Infield (2018), only consider the relationship between wind power and wind speed, and do not include additional variables that may affect these dynamic relationships and therefore affect the power curve accuracy and thus the effectiveness of power curve-based fault detection techniques. Recently, rotor speed and blade angle pitch (Singh, 2013) based data-driven models were used for turbine performance/condition monitoring, though the impact of these operational variables on power curve based data-driven models were not assessed. The study reported here incorporates these operational variables into power curve-based, data-driven machine learning models and analyses their impact on model accuracy and uncertainty.

SCADA data description, preparation and filtering (including air density corrections)
The SCADA system records more than 100 different signals, ranging from the timestamp, wind speed, wind direction, yaw angle, pitch angle, active power, reactive power, generator current, generator speed, gearbox temperature, generator winding temperature and ambient temperature. A SCADA dataset contains extensive real-time information that can be effectively applied to condition/ performance monitoring activities (e.g., failure or anomaly detection), most frequently 10-minute's average values.
In this research, data from a number of 2.3 MW, pitchregulated, variable speed onshore WTs located in Scotland, UK, have been obtained for model construction and validation purposes. Table 1 describes the characteristics of the six-month SCADA dataset, which starts with time stamp ''1 July 2012 00:00 AM'', ends at time stamp ''31 December 2012 23:50 PM'', and contains 30,628 measured values. Fig. 2 illustrates the measured power curve based on all 30,628 data points and suggests that a proportion of the SCADA data points recorded reflect measurement or recording errors. These erroneous data points can be due to sensor failures, missing data values or data collection faults and the data needs to be pre-processed to remove them as far as possible before further analysis. A practical yet straightforward SCADA data-based process where negative power values, mismatched timestamps, obscure (outlying) values and turbine power curtailment are identified is used to pre-process the data. However, it should be noted that the resulting SCADA data will   not be entirely free from errors, but the impact of residual errors is significantly reduced using these criteria. Following this process, the measured data points are reduced to 8256 after preprocessing. A small sample of the data used in this study is shown in Table 2, showing typical normal variation. The air density variations at wind farm site affects power curve accuracy, especially when a turbine is operating below rated power. Thus, it is usual to correct power curves based on air density before further analysis. The IEC 61400-12-1 standard acknowledges this and recommends that air density correction for pitch regulated wind turbines using the following equation: where B is the atmospheric pressure in mbar and T is the ambient temperature in • C, taken from the 10-minute averaged SCADA data values, as shown for example in Table 2. Thereafter, calculated air density values are used to obtain the corrected wind speed values using the following equation: Here, V C and V M are the corrected and measured wind speeds in m/sec.
The IEC prescribed air density corrections approach may not give the most accurate power curve in case of data-driven models, as pointed out by  where they proposed and carried out four possible air density correction approaches to GPbased power curves. In their study, they used extreme SCADA datasets of WTs located in extremely low-temperature (average monthly = −5.27 • C) and high-temperature (average monthly = 29.77 • C) regions both to comply with IEC standard guidelines that state that the air density correction shall be applied when the site density differs from the standard value (1.225 kg/m 3 ) by more than 0.05 kg/m 3 and to stress the influence of air density. The results suggest that the best option is to add air density directly into the GP model without pre-corrections. However, in this present paper, the SCADA dataset used is not from an extreme region and does not differ from the standard value by more than 0.05 kg/m 3 . Thus, it has been chosen to apply the standard IEC provisions for air density correction, as described above. Fig. 3 presents a pre-processed air density corrected power curve, which will be utilised in subsequent research and analysis.

Gaussian process methodology for power curve modelling
Due to the ease of operations and flexibility, there has been much activity concerning the application of GP to WT activities such as performance improvement and condition monitoring and they have been found to perform well (Pandit and Infield, 2018). A systematic and detailed explanation of the GP method can be found in Rasmussen and Williams (2006) so only a brief description is provided here as follows.
A GP model is a probability distribution over the function, which is defined by its mean function, , otherwise known as kernel. The f (x) is a GP distributed function are parameterised bym (x) and K ( x, x ′ ) using the following mathematical expression: The covariance (or kernel) function controls the correlation between all pairs of output values and therefore plays an important role in determining shape and smoothness of the function. There are various kinds of covariance functions available, whose performance varies depending upon the type of problem and so, choosing an appropriate kernel is important for developing a robust GP algorithm. Recently, authors of , carried out performance comparison of numerous kernels and concluded that the squared exponential is the most appropriate kernel for modelling power curves that is able to express the nonlinear relationship between wind speed and power output correctly. Hence, in this study, the squared exponential covariance function is chosen; it is mathematically defined as: A SCADA dataset of a wind turbine is not immune to noise, and therefore it is advisable to compensate for the effect of noise by adding a noise term into the covariance function, so Eq. (5) is rewritten as: where σ 2 f and l are defined as hyper-parameters. σ 2 f represents the signal variance, and l is a characteristic length scale which signifies how quickly the covariance decreases with the distance between successive data points. These hyper-parameters are interpretable and can be learned from data and optimised GP models to improve fitting accuracy, as described below.
The first step is to estimate mean and variance values from the SCADA training dataset to model the GP power curve. For instance, to predict the power y ′ for a new wind speed value x ′ for a training dataset, D = {(U i , P i ) , i = 1, . . . . . . , N}. These datasets are used to calculate the posterior distribution of P * for a given input U * and is defined as p(P * |U * , U tr , P tr ) in which {P * , U * } are the prospective power and wind speed values. The squared exponential covariance function depends on hyper-parameters which needs to be optimised before the posterior distribution of P * is calculated in order to ensure GP model accuracy. Maximisation of the log marginal likelihood is used to optimise hyperparameters (σ 2 f , σ 2 n and l), using the following equation (IEC Standard 61400-12-1, 2006), log (p (P tr |U tr )) = −0.5P T tr K −1 P tr − 0.5log (|K |) − 0.5nlog (2π ) (7) A quasi-Newton optimiser is used to optimise the hyperparameters concerning the likelihood in simple ML-II fashion,  (Rasmussen and Williams, 2006;Anon, 2019). After optimisation, estimation of the distribution of P * for a given U * is simple and straightforward.
The estimated distribution of P * , p (P * |U * , U tr , P tr ) follows a Gaussian distribution with mean and variance expressed by the following equations: is a covariance between test and training data points in the form of a column vector and, k * * = k (U * , U * ) is the auto-covariance function of the testing data points. The predicted mean of Eq. (8) is a continuous merger of the output P tr where linear weights are defined as k T * K −1 while the posterior variance of Eq. (9) is a function of k * , which is inversely proportional to the distance between test and training data points. It should be noted that, in this study, the 'OptimizeHyperparameters' function of matlab was used to calculate the hyperparameters; this minimises five-fold cross-validation loss automatically using Bayesian Optimisation, minimising the out-of-sample MSE as measured using cross-validation (Rasmussen and Williams, 2006). Rasmussen and Williams (2006) state the likelihood includes a trade-off between fit and model complexity, so overfitting tends to be less significant a problem in GP regression, but for complex system modelling, GP method are found to be complex as well as challenging. However, the main weakness is in dealing with large datasets in which every iteration of the optimisation requires the inversion of a N ×N matrix, where N is the number of data points. This leads to the asymptotic complexity issue called cubic inversion that is O (N) 3 and affects GP model accuracy. To deal with computational difficulties associated with large datasets, various methods have been proposed; see for example (Hartikainen and Särkkä, 2010). However, such methods need high processing power and have additional computational cost. Therefore, striking a balance between the number of data points and computation cost is important in order to construct a robust GP algorithm. For the above reasons, out of 8256 pre-processed data points, only 5000 data points are used for GP modelling out of which 3500 data points are used for training with the remaining 1500 data points used for testing purposes; a 70:30 ratio as illustrated in Table 3.
Using these datasets and above described methodology, a GP power curve has been calculated using and is shown in Fig. 4

Table 3
Data division for GP models training and testing.

Filtered data
Used data Training data Testing data 8256 5000 3500 1500 together with its intrinsic estimated 95% confidence intervals (CIs) demonstrating that the GP model is able to fit the power curve smoothly and accurately. It should be noted that the fitted GP model is constructed at this stage with only the density corrected wind speed used as input, i.e. without incorporating any operational variables.

Incorporating operational variables to gaussian process power curve models
As already described in the above section, the covariance function is the essence of the GP model and used to signify the similarity between two points. The general covariance matrix, K , gives the variance of each variable along the leading diagonal, and the off-diagonal elements measure the correlations between the different variables mathematically described as follows, where . K is of size n × n, where n is the number of input parameters considered, and it must be symmetric and positively semidefinite i.e. K ij = K ji. Due to multivariate nature of a general GP, n = number of predictors selected for GP model, while x represents the wind speed along with wind turbine performance parameters (rotor speed and rotor speed) to facilitate analysis of the effect of performance parameters on GP models. The operational variables are included individually along with wind speed to train the GP power curve model, as described in previous sections, and results are plotted in a 3D scatter plot to analyse the effect of incorporated variables on GP power curve. Uncertainty analysis plays a key role in assessing the performance of the models. GP models come with intrinsic confidence intervals (CIs) through its basis in Gaussian probability along with mean and variance estimates (see Eqs. (8) and (9)) and learn  the noise and smoothness parameters from the training data. CI values associated with each prediction, give vital information about uncertainty associated with the GP model and can be used to detect the outliers caused by wind turbine faults and failures as demonstrated by Pandit and Infield (2018). CIs are a good measure of the robustness and uncertainty of a model and can assist in understanding whether the inclusion of operational variables leads to any improvement in the GP model or not. The standard deviation is the square root of the variance of the predicted function (σ 2 (P * )) and used to estimate the CIs of the GP models using Eq. (11): CIs signify the pointwise mean plus and minus two times the standard deviation for given input value (corresponding to the 95% confidence region which represents the significance level of 0.05), for the prior and posterior respectively and will be used in future analysis for GP power curve uncertainty.

Incorporating blade pitch angle
The rotor power coefficient is a function of blade pitch angle (see Eq. (1)) which is used to regulate power production from a wind turbine. Averaging the angles of the three blades (in case of a three-bladed wind turbine) gives the value of mean pitch angle, and this is adjusted under normal operation to capture the maximum power below rated power and to limit power output above rated wind speed. Therefore, it affects the power production and shape of the power curve as well; for this reason, it is included as an extra input variable along with wind speed for GP power curve model training and testing purposes. Fig. 5 presents a 3D scattered plot of the estimated GP power curve, concluding that the inclusion of blade pitch leads to a slight improvement in accuracy. This is furthermore demonstrated in uncertainty analysis, where estimated CIs are plotted against corrected wind speed, as illustrated in Fig. 6, where results in improved uncertainty only occur above 13 m/s, compared to a GP power curve model without the inclusion of blade pitch angle. As already stated, the SCADA datasets are from pitch regulated wind  turbines which generally means that pitch adjustment is used only around and above rated wind speed range to limit power outputs in high winds, hence its impact on improving accuracy as well as uncertainty is not very significant, in fact its effect on the GP model can only be seen above a wind speed of 13 m/s, see Fig. 5 and Fig. 6. Therefore, blade pitch angle is different from other variables such as power and rotor speed and this difference is worth mentioning.

Incorporating rotor speed
Rotor speed is mostly used to establish the relationship between wind speed and rotor speed, called the rotor speed curve; and between rotor speed and power output, called the rotor power curve. These relationships are vital in detection of any anomaly/fault and further improving condition monitoring as demonstrated by Singh (2013). A GP-based power curve model trained with rotor speed included as an additional input along with wind speed. Then, the power curve is estimated, and compared with measured power curve as shown in Fig. 7. The estimated CI values (with and without the inclusion of rotor speed) are plotted as a function of corrected wind speed for uncertainty analysis purposes and are shown in Fig. 8. The inclusion of rotor speed into GP model reduces uncertainty significantly across the entire wind speed range as compared to the GP without rotor speed inclusion model, see Fig. 8. This also suggests that power curve-based GP fault detection for early fault detection should include rotor speed. It is worth noting that the results presented in this study indicate the importance of unsteady effects in turbine operation. If turbine tip speed ratio remained constant in the power-producing region, then there would be no added benefit from incorporating rotor speed in the model.

Comparative studies
In this section, a performance comparison of GP models is presented in which uncertainty analysis and performance error metrics are used to judge the impact of the additional performance parameters on GP model accuracy. Using GP power curve methodology and from Eq. (8), power values are estimated for a given wind speed. While doing so, performance parameters (rotor speed and blade pitch angle) are incorporated in order to analyse their impacts on estimated power values; the results are shown in Fig. 9. The impact of these additional parameters are clearly  shown, especially in the circled areas of Fig. 9 which suggest rotor speed makes the most significant improvement to GP model accuracy in estimating power values. However, the inclusion of blade pitch makes only a minor improvement, significantly less than rotor speed, which confirms the conclusions of previous sections.
Residuals are the difference between measured and predicted values and their distribution is informative in case of analysing the impact of operational variables on GP models. The frequency distributions of the calculated residuals of GP models incorporating rotor speed and blade pitch angle are shown in Fig. 10 together with a fitted Gaussian distribution and it is found that distribution of residuals for the GP with rotor speed included is close to Gaussian as compared to other GP models; suggest this model has relatively less bias. It should also be noted that in theory GP model residuals should be Gaussian distributed.
As described earlier, CIs are a useful measure of uncertainty and the precision of model estimates. Therefore, calculated CIs of GPs incorporating operational variables compared with estimated CIs values of a GP model without operational variables is plotted as a function of corrected wind speed, as shown in Fig. 11. Fig. 11 suggests that GP incorporating rotor speed improved the uncertainty significantly compared to the blade pitch angle. Furthermore, when rotor speed and blade mean angle are together incorporated into GP model, then its impact on improving uncertainty is highest, however, not far from the GP with only rotor speed, as shown in Fig. 11. This conclusion further demonstrated by the calculated numerical values of statistical error metrics as follows.

Quantifying GP model accuracy
To evaluate and compare the performance of the GP approaches in power curve estimation with the inclusion of operational variables, statistical error metrics are computed from both estimated and measured values. The most commonly used statistical error metrics include root mean square error (RMSE), determination coefficient (R 2 ) and mean absolute error (MAE) (Neyman, 1937). In this study, these three statistical error metrics are employed for evaluating the performance of the GP models. The mean absolute error used to measure the difference between the measured and predicted values, is estimated by the following equations: In terms of residuals, MAE = Where residuals are defined as : e = y − y ′ The root mean square error is the average of the squared error, is used to quantify the magnitude of the residuals and is mathematically defined as: In terms of residuals, this is: where y ′ corresponds to the GP estimated values for n different predictions and y to the measured values. The coefficient of determination (R 2 ) is calculated as the square of the correlation between predicted output and measured values (hence always in the range from 0 to 1 with values closer to 1 indicates better fitting of the model to the data) and is used to signify how close the data are to the fitted regression. It is defined as R 2 where SSE is the sum of squared errors, called 'residual sum of squares' and TSS is the total sum of squares, called total 'sum of squares'.
The calculated values of these performance error metrics for GP models are tabulated in Table 4. The calculated RMSE and MAE values for the GP with rotor speed reflect the amplitude errors and their small values indicate a strong relationship between measured and predicted values as compared to the GP with blade pitch angle. Furthermore, the GP with just rotor speed model gives R 2 = 0.9969; being so close to unity suggests a highly accurate model. However, the GP incorporating blade pitch angle and rotor speed model provides by a small margin the most accurate model. This statistical error metrics values are consistent with Fig. 11.

Conclusion
Wind turbines power curves are widely used in numerous WTs applications such as condition monitoring, wind power forecasting and prediction. The widely used IEC standard 'binning' method to calculate the power curve uses only the mean wind speed and air density as input variables. However, the power production of a WT is affected by other parameters, so taking these variables into power curve modelling can improve accuracy and thus improve the models' capabilities to detect early sign of failures.
In this study, a GP-based power curve model is proposed in which two vital operational variables, namely rotor speed and blade pitch angle are incorporated in order to evaluate their impact on GP model accuracy. Statistical error metrics, residuals analysis and uncertainty quantification are used to support the conclusions of the research. The GP power curve model with the inclusion of rotor speed yielded a significant improvement in both accuracy (Fig. 7) and uncertainty (Fig. 8). Furthermore, the inclusion of rotor speed makes distribution of residuals closer to Gaussian, than the others shown in Fig. 10. Furthermore, a comparative study suggests that the impact of rotor speed on GP power curve model accuracy and uncertainty is higher than for blade pitch angle, as shown by Fig. 11. This conclusion is further supported by calculated values of performance error metrics as summarised in Table 4.
Future work will extend this research to compare these results against other existing data-driven techniques and to apply the models in the development of effective and cost-effective condition monitoring techniques for wind turbines to support improved O&M strategies.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.