Bootstrapped Ensemble of Artiﬁcial Neural Networks Technique for Quantifying Uncertainty in Prediction of Wind Energy Production

: The accurate prediction of wind energy production is crucial for an affordable and reliable power supply to consumers. Prediction models are used as decision-aid tools for electric grid operators to dynamically balance the energy production provided by a pool of diverse sources in the energy mix. However, different sources of uncertainty affect the predictions, providing the decision-makers with non-accurate and possibly misleading information for grid operation. In this regard, this work aims to quantify the possible sources of uncertainty that affect the predictions of wind energy production provided by an ensemble of Artiﬁcial Neural Network (ANN) models. The proposed Bootstrap (BS) technique for uncertainty quantiﬁcation relies on estimating Prediction Intervals (PIs) for a predeﬁned conﬁdence level. The capability of the proposed BS technique is veriﬁed, considering a 34 MW wind plant located in Italy. The obtained results show that the BS technique provides a more satisfactory quantiﬁcation of the uncertainty of wind energy predictions than that of a technique adopted by the wind plant owner and the Mean-Variance Estimation (MVE) technique of literature. The PIs obtained by the BS technique are also analyzed in terms of different weather conditions experienced by the wind plant and time horizons of prediction.


Introduction
The contribution of wind energy to the electricity production portfolio is increasing compared to other productions with energy sources, such as nuclear, coal, hydroelectric, oil and gas, and biomass plants [1,2]. Wind energy production capacity grew by more than 18% (111 GW) in 2020, for an overall capacity of 733 GW [3]. In some EU countries, wind has become the mainstream source of electricity production, e.g., in Denmark, where it constitutes 48% of the country's total electricity consumption in 2020 followed by Ireland (38%), Portugal (25%), Germany (27%), and the UK (27%) [4].
Different prediction models have been developed and used to estimate energy productions in wind plants. Generally, we can distinguish model-based and data-driven approaches [12][13][14]. Model-based approaches employ physics models that utilize windforecasting data for wind energy predictions [8,[15][16][17]. Data-driven approaches do not use explicit physical models. Instead, they rely exclusively on wind data to build (black-box) obtained by other benchmarks for different wind farms located in Taiwan. Quan et al. [54] proposed a Particle Swarm Optimization (PSO)-based LUBE method for short-term load and wind power prediction and uncertainty quantification. The quality of the constructed PIs was superior to other PIs obtained by other benchmarks for the load and wind power production prediction in a short time.
In this work, the BS technique is considered. It has already been applied and shown to be effective for quantifying uncertainty in different industrial applications [41,47,49,57,58]. The PIs are obtained for a predefined confidence level α% [41], i.e., upper and lower bounds of the prediction that bracket the true/actual energy production value with probability α% [59].
Thus, in this work, an ensemble of ANN models is developed for energy production prediction, and the BS is used to quantify the uncertainties that affect the predictions.
BS is applied to a real case study of a wind plant located in the south of Italy with a capacity of 34 MW. For comparison, the technique (hereafter called the Quantile technique) adopted by the wind plant owner for computing the 10th and 90th percentiles of the predictions and the MVE technique of literature are considered. The application results show that BS is superior and more informative for the electric grid decision-makers than Quantile and MVE techniques.
The PIs obtained by the BS are also analyzed in terms of (i) the different operational conditions experienced by the wind plant (i.e., the weather conditions) and (ii) the time horizons (i.e., delays) of the predictions.
With respect to (i), the different weather conditions experienced by the wind plant are identified by resorting to Principal Component Analysis (PCA) [60,61] and Fuzzy C-Means (FCM) [62,63]. The former is used to reduce the dimensionality of the input weather-forecasting data, keeping the relevant information content. In practice, the F weather-forecasting quantities are projected into the space of F* < F principal components, which allow for describing a large fraction of the data variability using few features. The FCM algorithm receives, as an input, the identified F* principal components of the historical dataset and provides output clusters made by data characterized by similar weather conditions. The optimal number of clusters is identified by using the Davies-Bouldin (DB) validity index [64].
Once the different weather conditions are identified, the quality of the PIs obtained by the BS is quantified for each weather condition in terms of PI widths and PI coverage values.
With respect to (ii), the quality of the PIs obtained by the BS is evaluated for each day-ahead production prediction (i.e., until four days ahead) in terms of the Root Mean Square Error (RMSE) used for assessing the accuracy of the production predictions, and the widths of the corresponding PIs at each prediction hour.
The analysis shows that the PIs significantly vary due to the effects of the weather conditions and the time horizon of the predictions.
Thus, the significant contributions of this work are: • The quantification of the uncertainty affecting the wind energy production predictions provided by an ensemble of ANNs models employing BS technique; • The comparisons with the Quantile (adopted by the plant owner) and MVE (from the literature) techniques used for the uncertainty quantification; • The analysis of the BS's PIs in terms of various influencing factors, such as the different weather conditions experienced by the wind plant and the time horizons of the predictions.
The remainder of this paper is organized as follows. In Section 2, the ensemble of ANN models for predicting wind energy production is presented. In Section 3, the BS technique for constructing PIs is described. In Section 4, the results of applying the BS technique to a real case study of a wind plant are presented and compared with those obtained by the Quantile approach of the wind plant owners and the MVE technique in the literature. In Section 5, the PIs obtained by the BS technique are analyzed with respect to the different operational conditions experienced by the wind plant and the time horizons of the predictions. Finally, some conclusions are drawn in Section 6.

Ensemble Approach for the Prediction of Energy Production in Wind Plants
In this Section, the ensemble approach for predicting energy production is described. Ensemble approaches have been applied in various application fields to enhance the accuracy of the predictions and quantify their uncertainty [34,35]. The basic idea is that the individual models of the ensemble can complement each other by leveraging their strengths and overcoming their drawbacks: thus, the aggregation of their outcomes can boost the performance of the models [35,[65][66][67].
A typical ensemble of models for energy production prediction is shown in Figure 1. A training dataset X train , which comprises the weather-forecasting data (WF) and the associated energy productions ( → P), is utilized for building the individual models of the ensemble (hereafter called the base models). Once constructed, the ensemble provides the prediction of the energy production (P ensemble ) for any new coming test pattern, whose input weather-forecasting data → WF test j at time t j is provided by a weather-forecast model. For this, the ensemble performs two steps [35,66]: (1) energy production predictions by the individual prediction models, and (2) aggregation of the energy production predictions.

Ensemble Approach for the Prediction of Energy Production in Wind Plants
In this Section, the ensemble approach for predicting energy production is described. Ensemble approaches have been applied in various application fields to enhance the accuracy of the predictions and quantify their uncertainty [34,35]. The basic idea is that the individual models of the ensemble can complement each other by leveraging their strengths and overcoming their drawbacks: thus, the aggregation of their outcomes can boost the performance of the models [35,[65][66][67].
A typical ensemble of models for energy production prediction is shown in Figure 1. A training dataset , which comprises the weather-forecasting data ( ) and the associated energy productions ( ⃗ ), is utilized for building the individual models of the ensemble (hereafter called the base models). Once constructed, the ensemble provides the prediction of the energy production (̂) for any new coming test pattern, whose input weather-forecasting data ⃗⃗⃗⃗⃗⃗⃗ at time is provided by a weather-forecast model. For this, the ensemble performs two steps [35,66]: (1) energy production predictions by the individual prediction models, and (2) aggregation of the energy production predictions. Correspondingly, the two key components for constructing an effective ensemble approach are [35,66,68]: (1) a strategy for obtaining diversity among the base models and (2) a strategy for aggregating the outcomes of the base models. With respect to (1), different strategies have been developed for injecting the diversity among the base models of the ensemble, e.g., adopting different predictive modeling techniques, adopting the same prediction model type but with different parameter settings, and training each model with different training datasets, by resorting to techniques such as Bootstrapping AGGregatING (BAGGING) [69,70], Boosting [71], and Adaboost [72]. The reader interested in more details on the techniques used for generating diversity in the base models can refer to [69].
In this work, we use BAGGING to create the diverse base models of the ensemble. The basic idea of BAGGING is to train the base models with different training datasets generated by bootstrap [73]: the different versions of the training datasets are created by randomly sampling from the original training dataset with replacement. Correspondingly, the two key components for constructing an effective ensemble approach are [35,66,68]: (1) a strategy for obtaining diversity among the H base models and (2) a strategy for aggregating the H outcomes of the base models.
With respect to (1), different strategies have been developed for injecting the diversity among the base models of the ensemble, e.g., adopting different predictive modeling techniques, adopting the same prediction model type but with different parameter settings, and training each model with different training datasets, by resorting to techniques such as Bootstrapping AGGregatING (BAGGING) [69,70], Boosting [71], and Adaboost [72]. The reader interested in more details on the techniques used for generating diversity in the base models can refer to [69].
In this work, we use BAGGING to create the H diverse base models of the ensemble. The basic idea of BAGGING is to train the base models with different training datasets generated by bootstrap [73]: the different versions of the training datasets are created by randomly sampling from the original training dataset X train with replacement.
Artificial Neural Network (ANN) is employed as a base model for the prediction of energy production. The motivation for using an ensemble of ANNs is the fact that this model is currently used by some energy production companies for wind energy production prediction, and it has been shown capable of providing more accurate and robust predictions than individual models [38].
With respect to (2), different strategies have been proposed for an effective aggregation of the base models' outcomes into a final aggregated one, e.g., statistics-based (Simple Average (SA) and Simple Median (SM)), and model performance-based (Globally Weighted Average (GWA) and Locally Weighted Average (LWA) (Local Fusion (LF))) [35,74,75].
In practice, the aggregation of the base model outcomes entails (1) assuming a weight w h for the energy production predictionP h obtained by each base model h, h = 1, . . . , H, and (2) aggregating the H outcomes as a weighted average: SA weights all the outcomes of the base models with equal weights, i.e., w h = 1 H . SM takes only the center value of the H base model outcomes distribution, i.e., it assumes that the weights are all equal to 0 except for that corresponding to the median of the H base model outcomes. GWA and LWA consider weights for the base models that are inversely proportional to their prediction performance computed on a validation dataset and neighboring validation patterns close to the test pattern under study, respectively.

PIs for Uncertainty Quantification of Wind Production Prediction
A PI with confidence level α% is defined as an interval P lower j ,P upper j , such that the probability that the true/actual energy production value, P j , of the test pattern at the time t j falls within the interval is equal to α% [41,59]: For evaluating the estimated PIs, two indicators are considered: (i) the coverage, i.e., the fraction of the true/actual energy productions which actually fall within the constructed PIs, and (ii) the PIs width. A PI with confidence α% should have coverage of at least α%, with a width that is as small as possible [41,59,76].

Bootstrap (BS) Technique
In wind energy predictions, the prediction error variance, σ 2 ε , can be decomposed into three terms corresponding to the following three sources of uncertainty: • σ 2 WF is the variance caused by the model input uncertainty, i.e., the weather-forecast errors (source of uncertainty 1); • σ 2 PR is the variance caused by the stochastic variability of the physical process (source of uncertainty 2); • σ 2 MO is the variance caused by the ensemble model error, e.g., due to random initialization of the ANNs parameters or to the different datasets used for training the ANNs (source of uncertainty 3).
The prediction error variance, σ 2 ε , can be decomposed into the three contributions by: The flowchart of the BS technique for the estimation of the unknown σ 2 ε , and the associated PIs, is sketched in Figure 2. There are three steps: Step 1: Building the BS training dataset. Let us assume that we have available a dataset of weather-forecasting data and their associated energy productions, X all = WF all , This dataset is portioned into two datasets: a training dataset X train = WF train , valid caused by the ensemble model uncertainty, can then be estimated using X valid : Step 2: Constructing the BS PIs of the test pattern. The BS training dataset  formed by the weather-forecasting data of the validation dataset, WF valid , and the squared prediction errors of the ensemble on the energy productions of the validation dataset, j . Finally, the PI of the test pattern at a time t j with a confidence level equal to α% is [55,57]: whereP test j is the energy production predicted by the ANN ensemble for the test pattern at a time t j and C α do f is the (1 − α)/2 quantile of Student t-distribution with a number of degrees of freedom equal to the number of ensemble models H. ] for providing estimates of the energy productions, ⃗ , whose true/actual productions ⃗ are already known. The variance 2 caused by the ensemble model uncertainty, can then be estimated using : Step Finally, the PI of the test pattern at a time with a confidence level equal to % is [55,57]: where ̂ is the energy production predicted by the ANN ensemble for the test pattern at a time and is the (1 − )/2 quantile of Student t-distribution with a number of degrees of freedom equal to the number of ensemble models .

Case Study
In this section, the ensemble approach of Section 2 is applied to the estimation of the uncertainty affecting the prediction of wind energy productions based on available weather-forecasting data and corresponding known energy productions of a wind plant located in the south of Italy [

Case Study
In this section, the ensemble approach of Section 2 is applied to the estimation of the uncertainty affecting the prediction of wind energy productions based on available weather-forecasting data and corresponding known energy productions of a wind plant located in the south of Italy [38] (Section 4.1). The quantification of the three sources of uncertainty (namely, uncertainty due to the model input (weather forecasts), uncertainty due to the inherent variability (stochasticity) of the physical process, and uncertainty due to the model error) is carried out by the BS technique, described in Section 3.1, and the results are compared with those obtained by the technique adopted by the plant owner for the estimation of PIs and the Mean-Variance Estimation (MVE) technique of literature (Section 4.2).

Data Description and Ensemble Model Development
In this Section, the dataset of real weather-forecasting data, WF, and corresponding energy productions, → P, of a wind plant with 34 MW capacity is described. The dataset has been collected every three-hours over three years (from 2011 to 2013) with a forecast horizon ∆t = 96 h (four-day ahead). In other words, at a given time t, the weatherforecasting data of the following ∆t = 96 h are available, with a datum every 3 h, i.e., at time t, t + 3 h, t + 6 h, . . . , t + 96 h [38].
Engineering and expert judgment have been used to select a set of F = 19 features (whose detailed characteristics cannot be revealed, due to confidentiality reasons), e.g., wind speeds (in meters/second), horizontal (u) and vertical (v) wind components, hour which the weather forecasting is referred to, temperature, etc., for building the ensemble and predicting the energy productions. Note that for confidentiality reasons, throughout the paper, the values of the wind speeds and energy productions reported in Figures and Table are given on an arbitrary scale. Figure 3 presents the one-day ahead wind speed forecasts ( Figure 3a) and corresponding energy productions (Figure 3b) of the year 2013. Figure 3 shows the large variability in the plant's wind speed and the related large variability of the energy productions. Note that the wind speed sign (i.e., positive or negative) refers to the wind direction. For example, the negative wind speed values of the horizontal component (u) indicate that the direction of the wind is from west to east, whereas the negative wind speed values of the vertical component (v) indicate that the direction of the wind is from north to south.  The data are appended in the matrix X all , where rows and columns represent the forecasting patterns and the physical quantities of the weather forecasts with corresponding energy productions, respectively. The 2011-2012 data are divided randomly into X train training dataset (a fraction of 70% with N train patterns) and X valid validation dataset (remaining fraction of 30% with N valid patterns) to build the individual models and develop the ensemble, respectively. The 2013 data are used as a test dataset X test .
In this work, an ensemble composed of H = 100 ANNs models has been built. Each ANN is characterized by an architecture with four layers (one input, two hidden, and one output) and 9 × 7 hidden neurons, following a trials-and-errors procedure. Figure 4 shows two examples of energy production predictions (squares) and the corresponding true values (circles) for two different days of the year 2013. It can be seen that the estimated productions are reasonably close to the actual production values, although in some cases, the prediction error is not negligible (e.g., t = 6 h). In this work, an ensemble composed of = 100 ANNs models has been built. Each ANN is characterized by an architecture with four layers (one input, two hidden, and one output) and 9 × 7 hidden neurons, following a trials-and-errors procedure. Figure 4 shows two examples of energy production predictions (squares) and the corresponding true values (circles) for two different days of the year 2013. It can be seen that the estimated productions are reasonably close to the actual production values, although in some cases, the prediction error is not negligible (e.g., = 6 h).

Application Results of the BS PIs Estimation Technique
In this Section, the PIs obtained by the BS technique on the test data of the year 2013, , are presented and compared with those obtained by two other PIs estimation techniques: (1) the technique adopted by the wind plant owner (hereafter called the Quantile technique) and (2) the Mean-Variance Estimation (MVE) technique from the literature.
Briefly, the basic idea of the Quantile technique is to consider the quantiles of the predicted energy productions obtained by the ANN models of the ensemble of Section 2 at each time . The PIs obtained by the Quantile technique are made of the 10th and 90th percentiles (lower and upper bounds, respectively) of the energy production predictions obtained by the = 100 models of the ensemble, for a target confidence level = 80%.
With respect to the MVE technique, its basic idea is to assume that the prediction error obtained by the ensemble, i.e., = −̂, is an uncertain variable distributed according to a Gaussian distribution function whose variance 2 has to be estimated by using a dedicated ANN, adequately developed with a procedure similar to that carried out for the BS technique [55]. The dependence of this variance on the weather-forecasting data's input patterns is the fundamental assumption of the MVE (refer to Appendix A for more details on the MVE technique for PIs estimation [55]).
For each of the BS and the MVE techniques, a dedicated feedforward ANN, with an architecture of three layers (input, hidden, and output) and seven hidden neurons, is developed to estimate the prediction error variance with the BS and MVE techniques.

Application Results of the BS PIs Estimation Technique
In this Section, the PIs obtained by the BS technique on the test data of the year 2013, X test , are presented and compared with those obtained by two other PIs estimation techniques: (1) the technique adopted by the wind plant owner (hereafter called the Quantile technique) and (2) the Mean-Variance Estimation (MVE) technique from the literature.
Briefly, the basic idea of the Quantile technique is to consider the quantiles of the predicted energy productions obtained by the H ANN models of the ensemble of Section 2 at each time t j . The PIs obtained by the Quantile technique are made of the 10th and 90th percentiles (lower and upper bounds, respectively) of the energy production predictions obtained by the H = 100 models of the ensemble, for a target confidence level α = 80%.
With respect to the MVE technique, its basic idea is to assume that the prediction error obtained by the ensemble, i.e., ε = P −P, is an uncertain variable distributed according to a Gaussian distribution function whose variance σ 2 ε has to be estimated by using a dedicated ANN, adequately developed with a procedure similar to that carried out for the BS technique [55]. The dependence of this variance on the weather-forecasting data's input patterns is the fundamental assumption of the MVE (refer to Appendix A for more details on the MVE technique for PIs estimation [55]). For each of the BS and the MVE techniques, a dedicated feedforward ANN, with an architecture of three layers (input, hidden, and output) and seven hidden neurons, is developed to estimate the prediction error variance with the BS and MVE techniques. For illustration purposes, the energy production predictions (squares) and the true production values (circles) together with the PIs obtained by the BS (shaded area), Quantile (triangles), and MVE (diamonds) techniques for two different days of the year 2013 test data are shown in Figure 5. It can be seen that: For illustration purposes, the energy production predictions (squares) and the true production values (circles) together with the PIs obtained by the BS (shaded area), Quantile (triangles), and MVE (diamonds) techniques for two different days of the year 2013 test data are shown in Figure 5. It can be seen that:


The PIs obtained by the Quantile technique (triangles) are narrow but with very low coverage values, i.e., the actuals productions fall outside the estimated PIs;  The PIs obtained by the MVE technique (diamonds) are wider than the PIs of the Quantile technique, and, consequently, have larger PI coverage probabilities;  The PIs obtained by the BS technique (shaded area) are wider than the Quantile technique (triangles) and slightly wider than the MVE technique (diamonds), and, consequently, the BS PIs have larger PI coverage probabilities, i.e., the true productions fall inside the estimated BS PIs.

Factors Influencing the Estimated BS PIs
In this section, the PIs obtained by the BS technique are analyzed in terms of (i) the different weather conditions that influence the wind energy productions (Section 5.1) and

Factors Influencing the Estimated BS PIs
In this section, the PIs obtained by the BS technique are analyzed in terms of (i) the different weather conditions that influence the wind energy productions (Section 5.1) and (ii) the time horizons (i.e., delays) of the predictions (Section 5.2).

Influence of the Weather Conditions
The wind plant under study is, indeed, affected by a very large variability in the weather conditions (Figure 3). In practice, one might be interested in (1) identifying the different weather conditions that can be experienced by the plant, e.g., low, medium, and high wind speed values, based on the available weather-forecasting data, and consequently, (2) investigating their influence on the estimated BS PIs.
With respect to (1), the overall dataset, X all , has been further analyzed as follows: • The dataset is high dimensional, i.e., it comprises F = 19 physical quantities of the weather forecasts and, therefore, it has been transformed into F * fewer dimensions by resorting to Principal Component Analysis (PCA) [60,61].
Two principal components, F * = 2, have been selected as representative of the weather forecast data of the wind plant under study. Figure 6a shows that the selected PCs explain 97% of the original weather forecast data, whereas Figure 6b shows the overall dataset in the space of the identified principal components. The two selected principal components can describe the dataset with reasonable accuracy: indeed, Figure 7a shows that the two components are capable of reconstructing the original weather forecast data (for the sake of clarity, the first 500 h are plotted) with low reconstruction errors, i.e., with residuals close to 0 (Figure 7b).

•
In the space of the identified principal components, the dataset has been partitioned into S dissimilar groups (whose number is "a priori" unknown), such that data belonging to the same group are very similar to each other and dissimilar to those of the other groups. The S groups can be interpreted as different operating conditions of the wind plant that can influence the wind energy production.
To this aim, the data shown in Figure 6b are clustered by the unsupervised Fuzzy C-Means (FCM) algorithm [62,63]. For identifying the optimum number of the groups C opt , single clustering validity index (e.g., Silhouette, Davies-Bouldin (DB), etc.) or a combination of different validity indices can be used [77]. In this work, Davies-Bouldin (DB) validity criterion has been considered for clustering the groups of the dataset of Figure 6b. The Davies-Bouldin (DB) criterion is based on the ratio of within-group and between-group distances: the optimal partition, which gives optimal separation and compactness of the obtained groups, has the smallest DB index value [64].   With respect to (2), once the energy production predictions of the 2013 test data are obtained by the ensemble approach of Section 2, the corresponding PIs are estimated by the BS technique of Section 3.1 for the quantification of the uncertainties that affect the predictions. The estimated PIs are evaluated in terms of the three operating conditions of Figure 8b. Figure 9 shows the average PI widths (Figure 9a) and the average PI coverage values (Figure 9b) of the data of the three weather forecast groups. One can easily recognize that the larger the variability of the wind speeds (group 1 and group 3), the larger the prediction error, and, coherently, the larger the width of the PI (the values are indicated on Figure), but the smaller the PI coverage probability, and vice versa.
With respect to (2), once the energy production predictions of the 2013 test data are obtained by the ensemble approach of Section 2, the corresponding PIs are estimated by the BS technique of Section 3.1 for the quantification of the uncertainties that affect the predictions. The estimated PIs are evaluated in terms of the three operating conditions of Figure 8b. Figure 9 shows the average PI widths (Figure 9a) and the average PI coverage values (Figure 9b) of the data of the three weather forecast groups. One can easily recognize that the larger the variability of the wind speeds (group 1 and group 3), the larger the prediction error, and, coherently, the larger the width of the PI (the values are indicated on Figure), but the smaller the PI coverage probability, and vice versa. This can be explained by the fact that the larger the variability of the weather conditions (group 1 and group 3), the larger the wind energy production and, hence, the larger the uncertainty in the energy production prediction, as shown in Figure 9. For clarification purposes, Figure 9c shows examples of the estimated PIs of few data points of the three weather conditions.

Influence of the Time Horizon
The ensemble approach of Section 2 trained on the 2011-2012 data is used for the predictions of the energy productions of the 2013 test data for a time horizon (hereafter called delays) of four days, ∆ = 93 h; namely, delays 1-4 correspond to the predictions in the time intervals [0 − 21], [24 − 45], [48 − 69] and [72 − 93] , respectively. Figure 10a shows the average Root Mean Square Error (RMSE) used for evaluating the accuracy of the production predictions on the overall test data, whereas Figure 10b shows the average widths of the corresponding PIs at each prediction hour. One can easily recognize that the larger the time horizon of the prediction, the larger the ensemble prediction error, and, coherently, the larger the PI's width. This can be explained by the fact that the larger the variability of the weather conditions (group 1 and group 3), the larger the wind energy production and, hence, the larger the uncertainty in the energy production prediction, as shown in Figure 9. For clarification purposes, Figure 9c shows examples of the estimated PIs of few data points of the three weather conditions.

Influence of the Time Horizon
The ensemble approach of Section 2 trained on the 2011-2012 data is used for the predictions of the energy productions of the 2013 test data for a time horizon (hereafter called delays) of four days, ∆t = 93 h; namely, delays 1-4 correspond to the predictions in the time intervals [0 − 21], [24 − 45], [48 − 69] and [72 − 93], respectively. Figure 10a shows the average Root Mean Square Error (RMSE) used for evaluating the accuracy of the production predictions on the overall test data, whereas Figure 10b shows the average widths of the corresponding PIs at each prediction hour. One can easily recognize that the larger the time horizon of the prediction, the larger the ensemble prediction error, and, coherently, the larger the PI's width.
For clarification purposes, Figure 11 shows an example of the energy production predictions of 4 days ahead (i.e., ∆t = 93 h) with the corresponding BS PIs estimates. One can easily recognize that, for large production values, the PIs are enlarged to accommodate the large uncertainty that affects the predictions. In contrast, for small production values, the PIs are shortened due to the small uncertainty that affects the predictions.
As the last remark, the decomposition of the sources of uncertainty that affect the energy production predictions has shown that ( Figure 12 For clarification purposes, Figure 11 shows an example of the energy production predictions of 4 days ahead (i.e., ∆ = 93 h) with the corresponding BS PIs estimates. One can easily recognize that, for large production values, the PIs are enlarged to accommodate the large uncertainty that affects the predictions. In contrast, for small production values, the PIs are shortened due to the small uncertainty that affects the predictions. As the last remark, the decomposition of the sources of uncertainty that affect the energy production predictions has shown that ( Figure 12):


As expected, process and measure errors (circle) are increasing with the time horizon of the prediction due to the weather forecast errors;  Process and measure errors are following the variability of the electricity production. This explains the variability in the widths of the obtained PIs;  Model error (diamond) is stable with respect to the time horizon of the prediction;  The overall error (triangle) is consequently increasing with the time horizon of the prediction and following the variability of the electricity productions. For clarification purposes, Figure 11 shows an example of the energy production predictions of 4 days ahead (i.e., ∆ = 93 h) with the corresponding BS PIs estimates. One can easily recognize that, for large production values, the PIs are enlarged to accommodate the large uncertainty that affects the predictions. In contrast, for small production values, the PIs are shortened due to the small uncertainty that affects the predictions. As the last remark, the decomposition of the sources of uncertainty that affect the energy production predictions has shown that ( Figure 12):


As expected, process and measure errors (circle) are increasing with the time horizon of the prediction due to the weather forecast errors;  Process and measure errors are following the variability of the electricity production. This explains the variability in the widths of the obtained PIs;  Model error (diamond) is stable with respect to the time horizon of the prediction;  The overall error (triangle) is consequently increasing with the time horizon of the prediction and following the variability of the electricity productions.

Conclusions
In this work, we have considered the problem of quantifying the uncertainty that affects the predictions of the energy productions of a wind plant. The uncertainty quantification is carried out by constructing Prediction Intervals (PIs) with a predefined confidence level (e.g., 0.8). To this aim, the Bootstrap (BS) technique has been applied, and its capabilities have been verified on a real case study of a wind plant located in Italy. The obtained PIs have been evaluated by considering two indicators: (1) the coverage, i.e., the

Conclusions
In this work, we have considered the problem of quantifying the uncertainty that affects the predictions of the energy productions of a wind plant. The uncertainty quantification is carried out by constructing Prediction Intervals (PIs) with a predefined confidence level (e.g., 0.8). To this aim, the Bootstrap (BS) technique has been applied, and its capabilities have been verified on a real case study of a wind plant located in Italy. The obtained PIs have been evaluated by considering two indicators: (1) the coverage, i.e., the fraction of the true/actual energy productions which actually fall within the PIs, and (2) the PI width. Results show that the BS technique is superior and more informative for the electric grid operators than a technique based on the use of the quantiles of the ensemble model predictions, which is currently used by the plant owner, and the Mean-Variance Estimation (MVE). In practice, only the proposed method is able to cover within the prediction intervals a fraction of the true production values larger than the predefined confidence interval, which confirms its capability of properly describing all uncertainty sources. The PIs obtained by the BS technique have been further analyzed in terms of the different weather conditions experienced by the wind plant and the time horizon of the predictions. Future work will be devoted to (1) optimizing the PIs in order to obtain a trade-off between PI coverage and PI width, which is satisfactory for the decision-maker; (2) developing a framework for the effective use of the PIs in the formulation of energy bidding strategies.
Author Contributions: S.A.-D., P.B. and E.Z. were responsible for the conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing-original draft, writing-review and editing, visualization, supervision, and project administration, and M.L. was responsible for the resources, validation, project administration, and writing-review and editing. All authors have read and agreed to the published version of the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following acronyms/notations are used in this manuscript:

Appendix A The MVE Estimation Technique for PIs Estimation
The flowchart of the Mean-Variance Estimation (MVE) technique for the estimation of the unknown σ 2 ε , and the associated PIs, is sketched in Figure 1. , on the validation dataset.
Step 2: Constructing the MVE PIs of the test pattern. With the MVE training dataset, a dedicated feedforward ANN is developed for providing, at time t j , an estimate of the variance, σ 2 ε test j , associated with a general test pattern of weather-forecasting data, → WF test j . To ensure a strictly positive variance estimate, an exponential activation function is used.
Thus, a dedicated feedforward ANN characterized by an architecture of three layers (input, hidden, and output) and seven hidden neurons are developed to estimate the prediction error variance with the MVE technique.
Finally, the PI with a confidence level equal to α% of the test pattern at the time t j is obtained as per Equation (A1) [55,57]: whereP test j is the energy production predicted by the ANN ensemble for the test pattern at a time t j and C α do f is the (1 − α)/2 quantile of Student t-distribution with a number of degrees of freedom equal to the number of ensemble models H.