Next Article in Journal
Esoteric Pull and Esoteric Push: Two Simple In-Place Streaming Schemes for the Lattice Boltzmann Method on GPUs
Next Article in Special Issue
Modelling the Thermal Effects on Structural Components of Composite Slabs under Fire Conditions
Previous Article in Journal
Treetop Detection in Mountainous Forests Using UAV Terrain Awareness Function
Previous Article in Special Issue
Optimal Control of a Passive Particle Advected by a Lamb–Oseen (Viscous) Vortex
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Cluster-Based Analogue Ensembles for Hindcasting with Multistations

by
Carlos Balsa
1,*,†,
Carlos Veiga Rodrigues
2,†,
Leonardo Araújo
3,† and
José Rufino
1,†
1
Research Centre in Digitalization and Intelligent Robotics (CeDRI), Instituto Politécnico de Bragança, 5300-253 Bragança, Portugal
2
Vestas Wind Systems A/S, Vestas Technology Centre Porto, 4465-671 Leça do Balio, Portugal
3
Universidade Tecnológica Federal do Paraná, Campus de Ponta Grossa, Ponta Grossa 84017-220, Brazil
*
Author to whom correspondence should be addressed.
These authors contributed equally to this work.
Submission received: 29 April 2022 / Revised: 20 May 2022 / Accepted: 27 May 2022 / Published: 2 June 2022

Abstract

:
The Analogue Ensemble (AnEn) method enables the reconstruction of meteorological observations or deterministic predictions for a certain variable and station by using data from the same station or from other nearby stations. However, depending on the dimension and granularity of the historical datasets used for the reconstruction, this method may be computationally very demanding even if parallelization is used. In this work, the classical AnEn method is modified so that analogues are determined using K-means clustering. The proposed combined approach allows the use of several predictors in a dependent or independent way. As a result of the flexibility and adaptability of this new approach, it is necessary to define several parameters and algorithmic options. The effects of the critical parameters and main options were tested on a large dataset from real-world meteorological stations. The results show that adequate monitoring and tuning of the new method allows for a considerable improvement of the computational performance of the reconstruction task while keeping the accuracy of the results. Compared to the classical AnEn method, the proposed variant is at least 15-times faster when processing is serial. Both approaches benefit from parallel processing, with the K-means variant also being always faster than the classic method under that execution regime (albeit its performance advantage diminishes as more CPU threads are used).

1. Introduction

Short-term weather predictions by correlation with similar states in the past (analogues) were originally established by Lorenz [1], who suggested that two atmospheric states that are initially very close to each other will remain somewhat similar in the future. This was introduced as an alternative to classical weather forecasting based on systems of equations underlying deterministic Numerical Weather Prediction (NWP) models.
For many years, however, Lorenz’s proposal was discarded because of limited historical data on past weather conditions (especially over wide geographical areas) and insufficient computing capacity to implement his approach. Two decades later, van den Dool [2] revisited analogue-based short-range weather forecasting and found it to be feasible and effective when applied to limited geographical areas.
Monache [3] showed the applicability of an analogue scheme (named AN) for post-processing numerical weather forecasts to reduce systematic and random errors. The basic idea is that if previous forecasts (analogues) exist that are similar to the current NWP forecast (predictor), it is possible to produce a AN forecast by averaging the observations corresponding to these previous forecasts. The analogue prediction is then compared with the NWP prediction to infer the prediction error and thus improve the NWP forecast.
Later, Monache refined the use of analogues to estimate the probability distribution of the future state of the atmosphere [4]. Instead of focusing on improving a single deterministic NWP prediction, the goal was to derive a Probability Density Function (PDF) of the range of possible future states (forecasts)—a more realistic approach due to imperfect initial conditions and model limitations that generate prediction errors. In the same paper, the term analogue ensemble (AnEn) was coined to refer to both the observations that correspond to past analogue predictions and the method used to select these observations.
Since its introduction, the AnEn method has gained traction in many contexts, e.g., renewable energy management—in particular, wind and solar energy forecasting [5,6]. It has also been combined with other approaches, such as artificial neural networks, e.g., for predicting the power generated by photovoltaic power plants [7].
Studies aiming at the development of AnEn techniques favour applications in downscaling and forecasting rather than hindcasting, yet all topics share affinities. In the area of downscaling, we highlight the study conducted by Rozoff and Alessandrini comparing AnEn with convolutional neural networks (CNNs) for reconstructing high-resolution 10 m winds over complex terrain, where AnEn was found to produce lower errors than CNNs [8].
In the field of forecasting, AnEn is mainly used for the post-processing of forecast results. In recent years, post-processing with AnEn has been applied to the prediction of a wide range of meteorological variables and energy. For example, the AnEn technique has been used to improve the accuracy of ozone and particulate matter forecasts [9]. Solar forecasting, which involves the prediction of irradiance and solar power forecasting, has also attracted considerable attention in the last decade [10].
Among all post-processing forecast applications, there is a consensus that AnEn benefits from increasingly large training datasets. However, the large amount of data that must be processed to determine the analogues often makes the computational cost prohibitive [11]. Thus, there is a great interest in improving the computational efficiency of the AnEn method. One of the major outcomes of this interest was the development of the Parallel Analogue Ensemble (PAnEn) library [12]. This library, which provides an efficient parallel implementation of the AnEn method and user-friendly interfaces in R and C++, has been successfully applied to a huge dataset related to photovoltaic power generation [13].
This study aims to improve the computational efficiency of the AnEn method, not with the focus on the parallelization of the original algorithm but by reducing the number of operations required to compute the analogues, which is the most demanding task. The various computational loops needed to determine the analogue ensembles were replaced by a single step that employs clustering through the K-means method.
This new approach was developed in the field of hindcasting. Since this field has many similarities with downscaling and forecasting, the proposed cluster-based variant of the AnEn method could be easily adapted to these fields. Furthermore, the algorithm of the new approach is presented in detail. The effects of its most important parameters are examined with regard to the reliability and computational efficiency of the method.
This work is a revised and extended version of [14]. It provides a more complete mathematical and algorithmic formulation of the classical and K-means-based AnEn method. A complete characterization of the dataset used is given. The influence of additional method parameters is also investigated. The computational performance results cover a broader scenario of parallel execution and additional performance metrics are presented.
In the following, Section 2 reviews the classical AnEn method in the context of hindcasting; Section 3 formally presents the different ways of selecting analogues for the different approaches investigated in this work; Section 4 presents the meteorological dataset and the error metrics used to validate the predictions; Section 5 provides the numerical results obtained for different values of the main parameters of the K-means-based AnEn method; Section 6 is devoted to the comparison of the computational performance of the proposed method with the classical one; Section 7 finally concludes with an outlook on future work.

2. Hindcasting with the AnEn Method

The AnEn method plays an important role in weather hindcasting. Classically, weather hindcasting involves applying a forecast model to a past starting point to validate the model by comparing its forecasts with available observations (reanalysis). If some observations are not available, their time-series can be complemented (reconstruction) by using the corresponding forecasts as a substitute for the missing observations. However, the combination of hindcasting with the AnEn method also makes it possible to reconstruct meteorological data. For example, a variable at a meteorological station can be reconstructed based on the data of correlated variables from the same station or/and from other nearby stations.
The last option, known as hindcasting with multistations, was explored in previous work [15]. There, cosine similarity, normalization and K-means clustering were used as similarity metrics for analogue selection, as an alternative to the classical Monache metric. The coupling of the AnEn method with K-means clustering proved to be promising. At the same time, the results suggested the need to consider other clustering methods and to conduct a parametric study on important parameters of the methods under consideration.
Such study was initiated in [16] in the context of reconstructing a single meteorological variable. The study restricted the choice of analogues to the classical Monache metric, the K-means and the C-means clustering methods and confirmed that K-means provides the best accuracy. Heuristics were also identified for determining the number of clusters, the number of analogues and the analogue time span, which minimized the prediction errors.
Recently, the combination of K-means clustering with the AnEn method was further investigated in [14]. The mathematical formulation of the resulting approach was introduced to emphasize the important options and parameters of the new method. A parametric study was conducted in the context of the hindcasting of extra meteorological variables from the same dataset used in [16]. The computational performance of the K-means-based and classical AnEn methods were compared for a limited parallel execution scenario.
The application of the AnEn method to the reconstruction of missing data, from time-series of real-world observations, is illustrated by Figure 1, for the scenario of a single predictor station. In the figure, the historical dataset is a full record of past observations of a certain meteorological variable collected at the predictor station, and the observation dataset is an incomplete record of the same or a correlated variable at the predicted station (this record is complete for a training period but incomplete or absent for a prediction period).
The reconstruction of a missing value (predicted value) of the predicted station, concerning a certain instant t in time, unfolds in three different steps.
In step ①, a certain number of analogues are selected from the historical dataset based on the similarity of the past observations to a predictor at instant t. Both the predictor and each analogue are vectors of 2 k + 1 elements sampled at successive intervals of Δ t time-step. Each element is the value of a meteorological variable and k is a positive integer that represents the width of each half-window (into the past and future) around the central instant of the time window ( k = 1 in the scenario of Figure 1). Comparing vectors instead of single values considers the evolutionary trend of the meteorological variable around the central instant of the time window, allowing for the selection of analogues to be based on short-range weather patterns instead of single isolated events.
In step ②, the analogue map onto observations in the training period of the predicted station. This mapping is done only for the central instant of the analogue time window, meaning that, for each analogue vector, only a single observational value is selected.
Finally, in step ③, the observations selected are used to estimate the missing predicted value. If this value is indeed available as real observational data (as assumed in this paper), it then becomes possible to assess the prediction/reconstruction error.
In this work, the methodology described above is combined with K-means clustering as an alternative way of defining the set of analogues. In the ensuing text, the following notation is used: H stands for the historical time-series from where analogues are defined, O represents the measurements/observations time-series for the feature to be predicted, and P is the outcome prediction. Whilst O and P can be viewed solely as function of time t, when using multistations the history H will be an aggregate of time-series of multiple predictor stations, in which case, it will be a function also of the station s—that is, H = H ( s , t ) .

3. Searching for Analogue Ensembles

This section presents the techniques that are relevant in the scope of this work to find analogues. It begins with the classical error-based approaches and ends with the proposed cluster-based technique (for which a complete formal description is provided).

3.1. Classical Error-Based Analogues

The classical techniques are based on the score obtained through given error metrics ( ϵ ). These metrics differ depending on the predictor stations being used in an independent or dependent way: in the first scenario, the choice of analogues for a predictor station is independent of the choice made for other predictor stations; in the second scenario, the analogues chosen for the different predictor stations must coincide in time.

3.1.1. Independent Analogues

With a single predictor station, a univariate similarity metric can be defined [3] as
ϵ ( t , τ ) = r = k k H ( t + r Δ t ) H ( τ + r Δ t ) 2
where H is a single historical time-series from where analogues are chosen, t is a time instant in the prediction period of H, τ is a time instant defining a possible analogue in the training period of H, and Δ t is the time-series time-step. Both t and τ are the central instants of time windows with 2 k + 1 consecutive instants that are Δ t apart in time. For each time window, there is a vector with 2 k + 1 consecutive records of a meteorological variable. The metric ϵ ( t , τ ) is a Euclidean distance between a possible analogue and the predictor, such that the best analogues are those that yield the lowest values (scores) of ϵ ( t , τ ) . This metric can be adapted to a scenario with multiple predictor stations at different locations, identified by an index s, each with its own time-series H ( s ) , turning into
ϵ I ( s , t , τ ) = r = k k H ( s , t + r Δ t ) H ( s , τ + r Δ t ) 2
where s = 1 , , N s , and N s is the total number of predictor stations. The metric ϵ I ( s , t , τ ) allows us to identify the best analogues for each station s, independently of those found for other stations; hence, this metric is designated as an independent score. This is equivalent to finding, for each predictor station s, the N a analogues with the lowest scores, i.e.,
a s , n = argmin τ ϵ I ( s , t , τ ) if τ { a s , 1 , , a s , n 1 } for s = 1 , , N s and n = 1 , , N a .
The prediction follows from the arithmetic mean of the observations in the time-series O, corresponding to the central instant of the time window of the selected analogues:
P ( t ) = 1 N s N a s = 1 N s n = 1 N a O ( a s , n ) .
This assumes the same number of best analogues being considered in each predictor station. However, it is also possible to consider, for each station, a different amount (c.f. [15,17]).
Algorithm 1 sums up the classic AnEn method with independent stations. The most demanding part from the computational point of view corresponds to the inner loops between lines 3 and 6; in them, the metric ϵ I is calculated m times per each predictor station, where m = N τ 2 k , and N τ is the total number of records in the training period of H, a period that may span several years (the number of subsets is m = N τ 2 k because each one of them has dimension 2 k + 1 and is constituted by the central register H ( τ ) plus the k previous and the k following records, therefore not being possible to form subsets for the first k and for the last k records of H); once these inner loops are repeated for each prediction i (outer loop between lines 2 and 12), the computational demand increases further, in direct proportion to the overall number of predictions, N p . However, this algorithm is easily parallelizable due to the data-independence of the iterations operations in its various loops.
Algorithm 1: Classic AnEn for independent stations.
1.
inputs: H, O, N p , N a , m and k
2.
for i = 1 , , N p
3.
for s = 1 , , N s
4.
  for j = 1 , , m
5.
    ϵ I ( s , t i , τ j ) r = k k H ( s , t i + r Δ t ) H ( s , τ j + r Δ t ) 2
6.
  endfor
7.
  for n = 1 , , N a
8.
    a s , n argmin τ ϵ I ( s , t i , τ ) if τ { a s , 1 , , a s , n 1 }
9.
  endfor
10.
endfor
11.
P ( t i ) 1 N s N a s = 1 N s n = 1 N a O ( a s , n )
12.
endfor

3.1.2. Dependent Analogues

When using several predictor stations, the analogues from different stations can be forced to overlap in time, meaning that there is a time dependency between the analogues. In this scenario, the score metric for the same instant τ in all time-series H ( s ) is given by
ϵ D ( t , τ ) = s = 1 N s ϵ I ( s , t , τ ) 2 .
Overall, there are now only N a best analogues to identify, and they are given by
a n = argmin τ ϵ D ( t , τ ) if τ { a 1 , , a n 1 } for n = 1 , , N a .
The best analogues thus correspond to historical vectors from different stations that coincide in time and, when considered together in Formula (5), ensure the lowest ϵ D values.
In turn, this translates into predictions based only on N a observations, given by
P ( t ) = 1 N a n = 1 N a O ( a n ) .
Algorithm 2 describes the classic AnEn method with dependent stations. Apparently simpler, this algorithm ends up involving similar computational effort as Algorithm 1 used with independent stations, as it also requires ϵ I to be computed m × N s times.
Algorithm 2: Classic AnEn for dependent stations.
1.
inputs: H, O, N p , N a , m and k
2.
for i = 1 , , N p
3.
for j = 1 , , m
4.
   ϵ D t i , τ j s = 1 N s ϵ I s , t i , τ j 2
5.
endfor
6.
for n = 1 , , N a
7.
   a n argmin τ j ϵ D t i , τ j if τ j a 1 , a 2 , , a n 1
8.
endfor
9.
P ( t i ) 1 N a n = 1 N a O ( a n )
10.
endfor

3.2. Cluster-Based Analogue Ensembles

The search for analogues in the training period of the historic time-series H may take considerable time due to (i) the need to go through every instant τ j in the training period and (ii) compare the vector of records centred in that instant with the predictor vector in order to compute the metrics ϵ I or ϵ D . In addition, parallelization (an alternative fast method to define the best analogues) was achieved by employing clustering techniques. Next, a description of this approach is provided, assuming the usage of K-means clusterisation.

3.2.1. Independent K-Means Analogues

The historic time-series H ( s ) of a station s may be broken into smaller overlapping vectors or subsets, of size 2 k + 1 , such that each subset j of the station s is given by
x s . j = { H ( s , ( j k ) Δ t ) , , H ( s , ( j + k ) Δ t ) } for s = 1 , , N s and j = k + 1 , , N τ k
where N τ is, as already stated, the dimension of the time-series of the training period. There are thus N τ 2 k subsets (vectors) x j per station. The set of these subsets for a station s is
X s = { x s . j } = { x s . ( k + 1 ) , , x s . ( N τ k ) } .
The clustering method is a function f that maps X s into a set of clusters c s . q , for a maximum number of clusters N c —that is, f : X s { c s . q } , with q = 1 , , N c , or identically
f : X s { c s . 1 , , c s . N c }
where each cluster c s . q will include a certain number of x s . j subsets that share an aggregation criteria (for instance, minimizing the coherence within each cluster as will be stated later). The aggregation of the x s . j subsets into a cluster will thus depend on the clustering algorithm employed and respective efficiency metric used. The number N c of clusters to be formed may be specified a priori, or estimated from X s , depending on the technique adopted.
After the application of the clustering algorithm, each cluster c s . q will have a centroid c ¯ s . q , corresponding to the mean of the cluster subsets (vectors). This centroid is given by
c ¯ s . q = x s . j c s . q x s . j x s . j c s . q 1 for q = 1 , , N c where c ¯ s . q c ¯ s . q ( τ ) τ [ k Δ t , k Δ t ]
meaning the centroid of a cluster is a vector and each element of that vector is given by the average of the corresponding elements in the vectors that belong to the cluster.
Each centroid vector c ¯ s . q acts as an individual analogue that may be compared against the historic value H ( s , t ) for a prediction time t, using a metric similar to Formula (2):
ϵ C ( s , t , c ¯ s . q ) = r = k k H ( s , t + r Δ t ) c ¯ s . q ( r Δ t ) 2 .
Having ranked all clusters of a station s by the Euclidean distance of their centroids to H ( s , t ) , it becomes possible to select the N a c best clusters. These will be the clusters c s . q whose centroids c ¯ s . q ensure the N a c lowest values of ϵ C ( s , t , c ¯ s . q ) :
c s , q = argmin c ¯ s . q ϵ C ( s , t , c ¯ s . q ) if c ¯ s . q { c s , 1 , , c s , q 1 } for s = 1 , , N s and q = 1 , , N a c .
Moreover, for each of the N a c clusters selected, one may consider all its members (vectors) as analogues—the approach adopted in this work—or only the N a best (the N a vectors closest to the cluster centroid, based on the clustering algorithm used).
Each subset (vector) x s . j of a cluster c s . j has a time correspondence to the observation time-series O that can be mapped into a matching subset of observations o s · j , given by
o s · j = O ( s , j Δ t ) for s = 1 , , N s and j = k + 1 , , N t k
It follows that each centroid c ¯ s . q will have an associated observation o ¯ s . q , which is the average of all the observations o s . j that matches the central time of the vectors x s . j c s . q :
o ¯ s . q = x s . j c s . q o s . j x s . j c s . q 1 for q = 1 , , N a c .
Considering the contribution of all predictor stations, the prediction is thus given by
P ( t ) = 1 N s N a c s = 1 N s q = 1 N a c o ¯ s . q .
In this work, a single analogue cluster is used per predictor station: c s . 1 , the cluster with the best score ϵ C ( s , t , c ¯ s . q ) . Thus, the prediction with N s predictor stations is simply
P ( t ) = 1 N s s = 1 N s o ¯ s . 1 .
The full sequence of the steps described above can be found in Algorithm 3. Clusterisation is performed only once (lines 2 to 4), for each of the historical datasets (one dataset per predictor station). Compared to Algorithm 1, the performance advantage of Algorithm 3 lies in the fact that the metric ϵ C ( s , t i , c ¯ s . q ) is computed (line 8) for a number of clusters N c that is usually much smaller than the number m of vectors for which the metric ϵ I ( s , t i , τ ) must be calculated in Algorithm 1 (line 5). Moreover, as in the classical algorithm, the main loops of the cluster-based algorithm may also be easily parallelized.
Algorithm 3: Cluster-based AnEn for independent stations.
1.
inputs: X s , H, O, N p , N c , N a and k
2.
for i = 1 , , N s
3.
{ c s . 1 , , c s . N c } f X s
4.
endfor
5.
for i = 1 , , N p
6.
for s = 1 , , N s
7.
  for q = 1 , , N c
8.
    ϵ C ( s , t i , c ¯ s . q ) r = k k H ( s , t i + r Δ t ) c ¯ s . q ( r Δ t ) 2
9.
  endfor
10.
  for q = 1 , , N a c
11.
    c s , q argmin c ¯ s . q ϵ C ( s , t , c ¯ s . q ) if c ¯ s . q { c s , 1 , , c s , q 1 }
12.
    o ¯ s . q x s . j c s . q o s . j x s . j c s . q 1
13.
  endfor
14.
endfor
15.
P ( t i ) 1 N s N a c s = 1 N s q = 1 N a c o ¯ s . q
16.
end

3.2.2. Dependent K-Means Analogues

The previous approach may be classified as independent in the sense that the historical dataset of each predictor station s is clusterised autonomously. As a result, the vectors of the best clusters of each station are not required (neither expected) to be perfectly aligned in time or even to overlap. It is however possible to apply clustering in a way that enforces some temporal correlation between the vectors of different stations, as next described.
Start by joining the vectors x s · j for the same j and different stations s. This produces a new vector x j with ( 2 k + 1 ) N s elements:
x j = s = 1 N s x s . j with j = k + 1 , N τ k .
The new vector x j may be viewed as a stripe made of slices x s · j . The set of all stripes is
X = { x j } = { x k + 1 , , x N τ k } .
Now, the clustering algorithm f will create N c clusters, bringing together correlated stripes (which implies their composing slices must also be correlated), into the same cluster:
f : X { c 1 , , c N c } .
Each cluster c q will have its own centroid c ¯ q . This centroid is still a vector of averages; however, with dependent stations, it has ( 2 k + 1 ) N s elements, being given by:
c ¯ q = x j c q x j x j c q 1 for q = 1 , , N c .
In turn, the predictor vectors H ( s , t ) are also joined, by the same order the vectors x s · j were joined; thus, there is now a single unified predictor vector, H ( t ) , with ( 2 k + 1 ) N s elements:
H ( t ) = s = 1 N s H ( s , t ) .
The predictor vector H ( t ) is compared against all N c centroids, using the metric:
ϵ C ( t , c ¯ q ) = r = k k H ( t + r Δ t ) c ¯ q ( r Δ t ) 2 .
Only the N a c best clusters/centroids are selected, being the ones that minimize ϵ C ( t , c ¯ q ) :
c q = argmin c ¯ q ϵ C ( t , c ¯ q ) if c ¯ q { c 1 , , c q 1 } for q = 1 , , N a c .
Once a set of clusters is chosen, the selection of the observations and the calculation of the prediction proceeds similarly to the way in which is performed for independent stations, except that now the clusters and the corresponding observations are not iterated by station:
o ¯ q = x j c q o j x j c q 1 for q = 1 , , N a c .
P ( t ) = 1 N a c q = 1 N a c o ¯ q .
Algorithm 4 resumes the main steps of the cluster-based AnEn method, considering dependent stations. Compared to Algorithm 3, clustering is done only once, irregardless of the number of predictor stations, albeit clusterisation must be preceded by the unification of their datasets. Thus, depending on the size and original organization of these datasets, such reorganization may offset the gains of a single clusterisation process (which will also have to deal with N s times more data).
There are also fewer clusters formed than in the independent stations scenario, and thus the error metric will need to be computed fewer times. The same applies to the number of average observations that need to be calculated once there are fewer centroids to consider. It remains to be seen, however, if the cluster-based dependent approach is faster than the independent variant, as in this work, these two approaches were only compared concerning their accuracy. Algorithm 4 should be, nevertheless, faster than its classic counter-part (Algorithm 2), due to the reasons already discussed that make clusterisation approaches inherently faster than classic ones.
Algorithm 4: Cluster-based AnEn for dependent stations.
1.
inputs: X , H, O, N p , N c , N a and k
2.
{ c 1 , , c N c } f X
3.
for i = 1 , , N p
4.
for q = 1 , , N c
5.
   ϵ C t i , c ¯ q r = k k H ( t i + r Δ t ) c ¯ q ( r Δ t ) 2
6.
endfor
7.
for q = 1 , , N a c
8.
   c q argmin c ¯ q ϵ C ( t , c ¯ q ) if c ¯ q { c 1 , , c q 1 }
9.
   o ¯ q x j c q o j x j c q 1
10.
endfor
11.
P ( t i ) 1 N a c q = 1 N a c o ¯ q
12.
endfor

3.2.3. K-Means Clusterisation

In the cluster-based AnEn method developed in this work, clustering is achieved using the K-means algorithm, whereby m = ( N τ 2 k ) N s historical data subsets x j I R 2 k + 1 are to be classified in N c clusters. The data is organized as lines in a matrix X I R m × ( 2 k + 1 ) . To describe the K-means method as proposed in [18], a partition of the subsets vectors x 1 , , x m in N c clusters is denoted as = c 1 , , c N c , where
c j = : x cluster j
defines the set of vectors in cluster j. The centroid (arithmetic mean), of the cluster j is then
c ¯ j = 1 n j c j x
where n j is the number of elements in cluster j. The sum of the squared distance, between the vectors and the j cluster centroid, known as coherence, is
q j = c j x c ¯ j 2 .
The closer the vectors are to the centroid, the smaller the value of the coherence q j . The quality of a clustering process can be measured as the overall coherence:
Q = j = 1 N c q j .
K-means is an optimization method. It searches for an optimal partition that minimizes Q ( ) . The problem of minimizing overall coherence is NP -hard [19], and consequently there is no guarantee that it will converge to the global optimum. The basic algorithm for K-means clustering is a two-step heuristic procedure. First, each vector is assigned to the group that is closest to it. Then, new centroids are calculated based on the assigned vectors. In the K-means version of Algorithm 5, adapted from [18], these steps are alternated until the changes in the overall coherence are smaller than a certain predefined tolerance.
Algorithm 5: K-means Algorithm.
  • Start with an initial random partitioning ( 0 ) and compute the corresponding centroid vectors c ¯ j ( 0 ) for j = 1 , , N c . Compute Q ( ( 0 ) ) . Set z = 1 .
  • For each x i find the closest centroid. If the closest centroid is c ¯ p z 1 assign x i to c p ( z ) .
  • Compute the centroids c ¯ j ( z ) for j = 1 , , N c of the new partitioning ( z ) .
  • If Q ( ( z ) ) Q ( ( z 1 ) ) < tolerance , stop; Else z = z + 1 and return to step 2.
As the initial partition is generated randomly, each execution of the K-means algorithm can lead to a different solution, all of them corresponding to quasi-optimal solutions that verify the convergence criterion.

4. Meteorological Dataset and Prediction Error Metrics

This section introduces the meteorological dataset used for the experimental evaluation and the error metrics used to asses the accuracy of the proposed hindcasting method.

4.1. Meteorological Dataset

All experiments presented in this work were conducted based on data from three meteorological stations located on the coast of the state of Virginia (USA): YKTV2, YKRV2 and DOMV2. Data for these and many other stations is freely available from the United States National Data Buoy Center [20]. The location of the selected stations is shown in Figure 2.
In a previous work [15] the correlation between the meteorological data of the different stations was analysed and the data series of station YKTV2 was reconstructed from one and two of the remaining stations, using different metrics in addition to the classical metric proposed by Monache [3] to find the analogues. It was then found that the weather at stations DOMV2 and YKRV2 followed a similar pattern although, due to their different locations, certain meteorological values differ.
However, since the alternative metrics look for a similar evolution of the weather during the time window and not for similar values as the classical metric, both stations are suitable for data reconstruction at station YKTV2. Therefore, these three stations were also used for the same purpose in this paper. Table 1 contains their coordinates and roles.
The data used in this study range from 2011 to 2019 and includes six different meteorological variables: Pressure (PRES), Air Temperature (ATMP), Mean Wind Speed (WSPD) and Gust Speed (GST). PRES is the pressure [hPa] measured at sea level, ATMP is the air temperature [°C], WSPD is the wind speed [m/s] averaged over an 8 min period for buoys and a 2 min period for land stations, and GST is the peak gust speed of 5 or 8 s gust speed [m/s] measured during the 8 min or 2 min period (for details see [20]). The time-series records are based on samples taken every 6 min ( Δ t = 360 s ). The time-series are not always complete. Data are missing for periods that can be short (hours) or long (months)—Table 2 shows the availability of data over the period studied.
Table 2 contains basic statistics of the relevant variables for each station, such as the minimum (Min), average (Mean) and maximum (Max) values, along with the number of missing records (# NA) and the percentage of available records (Availability) in the time-series. We verified that the amount of missing data did not exceed 3.3%. If the missing records did not exceed the size of the subsets x j , they were filled by interpolation from nearby values. If the range was longer, the corresponding records were ignored.
The 9 years of data was divided into two groups: (i) data for a training period, from the beginning of 2011 to the end of 2017; (ii) data for a prediction period, from the beginning of 2018 to the end of 2019. In the prediction period, data for station S 1 was hindcast from data of the stations S 2 and S 3 , hence acting as predictor datasets.
Having data for a period of 9 years with a time-step Δ t = 6 min implied a large number of data set vectors, which caused a considerable computational effort when applying the classical AnEn method. Therefore, the predictions in the forecast period were limited to the interval between 10 am and noon, based on the same intervals in the historical dataset. Theoretically, 10 samples per hour, 2 h per day, for 7 historical years would correspond to a total number of 10 × 2 × 365 × 7 = 51,100 samples and, accordingly, to the same number of possible analogue vectors (each centred in one of these samples); in the end, only a total of N τ = 43,238 records were used due to some missing data.

4.2. Prediction Error Metrics

As real data was available for the predicted period, it was possible to compare the predictions ( p i ) with the observed values ( o i ). Following Chai and Draxler [21] recommendations, multiple metrics were used to assess the model accuracy, starting with the Bias,
Bias = 1 N p i = 1 N p ( p i o i ) ,
where N p is the number of predictions, p i is a prediction, and o i is the corresponding truth value. The Bias measures the average error compared to the truth, however, allowing over- and under-predictions to cancel out. This metric can be interpreted as a rough approximation of the systematic error in the prediction. Thus, complementary to the Bias, the Root Mean-Squared Error (RMSE) is also used:
RMSE = 1 N p i = 1 N p ( p i o i ) 2 .
The RMSE is useful because the squared terms give a higher weight to higher errors. Thus, the RMSE will be higher if the model makes predictions which are far from the truth, even if these erroneous predictions are few in number.
The Standard Deviation of the Error (SDE) is computed from the Bias and RMSE metrics, corresponding to:
SDE = RMSE 2 Bias 2 .
The SDE represents the error due to variance; hence, it can be used as a rough approximation of the random error in the prediction.
Another metric, akin to the RMSE, is the Mean Absolute Error (MAE), being recommended by Chai and Draxler [21] as a complement to the former. The MAE is defined as
MAE = 1 N p i = 1 N p | p i o i |
thus, computing the average distance to the truth in absolute values. This is different from the Bias as the errors do not cancel out, whilst yielding a value smoother than the RMSE.
A low Bias and a high MAE indicates that the model is not really accurate but that its predictions are sometimes higher and sometimes lower than the truth. Thus, considering the MAE is necessary to really understand how the error in the prediction is distributed, as it also shows the systematic error but this time in terms of the absolute distance.
In addition to these individual metrics, this paper also uses a combined error metric (CE): C E = | B i a s | + R M S E + S D E + M A E , which is the sum of the absolute values of the individual error metrics (note that only Bias can be negative). When evaluating the impact of different values of critical parameters on the methods under study, the combined metric is a heuristic that allows the most favourable (in terms of error) parameter values to be easily identified: these values are those that lead to the absolute minimum of CE.
However, since many of the values of CE registered in the experiments are close to each other, a small tolerance of up to 10 2 is allowed; this means that values of CE within an excess of up to 10 2 from the absolute minimum are considered statistically identical, leading to several equivalent choices for the parameters of the methods under investigation.

5. Experimental Evaluation

In this section, the K-means AnEn method is evaluated and compared with the classical AnEn approach, in the context of a hindcasting problem with the previously presented meteorological dataset. All results were produced by implementations of both AnEn approaches in R [22]. For the K-means variant, the built-in function kmeans was used.

5.1. Variation of the Number of Clusters

The evaluation starts by analysing the impact of the variation of the number of clusters formed ( N c ) in the prediction errors. This number, and the amount of clusters effectively used ( N a c , with N a c N c ), are fundamental parameters of the cluster-based approach to the AnEn method investigated in this work: the number of clusters formed has an influence on the number of subsets (vectors) in each cluster (fewer clusters will have more subsets, and vice versa); correspondingly, this determines the number of observations used to derive the predictions and, ultimately, that number will affect the predictions accuracy.
The assessed scenario concerns the dependent stations variant, where a single cluster ( N a c = 1 ) is selected in each predictor station to generate the forecast. The results, shown in Figure 3, were obtained with k = 5 ; thus, for each meteorological variable there were m= N τ 2 = 43,228 subsets (vectors) x j in each predictor dataset from which clusters could be formed and each subset consisted of 11 values (covering a time frame of 66 min).
In Figure 3, a similarity in the behaviour of the errors with the variation of N c for the variables ATMP, GST and WSPD can be seen: the errors first decrease significantly with the increase of N c up to N c 50 ; then they decrease moderately up to N c 250 and stabilize thereafter. The variable PRES shows a different behaviour: the errors continue to decrease with the increasing of N c . This is possibly due to the fact that the variable PRES does not show large fluctuations over a short period of time, and therefore the smaller clusters describe its behaviour best due to the lower variance.
Determining the optimal value of N c a priori is difficult because each variable has different behaviour. Setting N c bellow 100 may result in a high error for each variable. For the variables ATMP, GST and WSPD, there is a small tendency for errors to increase as the value of N c increases, which would justify rejecting large values for N c . However, this is not true for the variable PRES, as the errors are inversely proportional to the N c value.
A possible heuristic to define a single common N c value is to set it as the square root of the total number of records in the training period ( N τ ) [16]. Considering the datasets used in this work, it would follow that N c = 43,238 208 . Another approach is to define it as the average of the N c values that ensure the lowest prediction errors for the different variables. In this case, the variable PRES is left aside (due to its singular behaviour) and, based on the values shown in Figure 3 for the other three variables, a “universal” value for N c would be ≈350; although this is not an optimal value for the variable PRES, it is a good compromise, since the error for this variable decreases only slightly after N c = 350 .

5.2. Variation of the Number of Analogues

Having determined the number of clusters to be formed ( N c ) and effectively used ( N a c ), there is still the possibility of refining the prediction by limiting the number of analogues N a (i.e., the number of x j subsets or vectors) provided by each selected cluster. Thus, it is possible to use all the subsets x j contained in a selected cluster or, alternatively, a smaller set ordered by affinity to the centroid of that cluster.
The effects of varying N a on forecast errors can be seen in Table 3 for the reconstruction of the four meteorological variables studied. These results were obtained keeping some parameters of the evaluation of the previous section (dependent stations, k = 5 , N a c = 1 ) but now with a limitation of the number of clusters formed to N c = 350 (according to the heuristic defined in the same section). In the table, the rows with the smallest verified errors for each of the meteorological variables are highlighted in bold. Moreover, the special value “∞” for N a means that all analogue candidates (all subsets/vectors) of a cluster were used.
The results in Table 3 show that, with the exception of the pressure variable (PRES), the lowest error values were obtained at high N a values. However, the differences between the smallest errors and those obtained with all subsets of the cluster were very small. This result indicates that for the variables ATMP, WSPD and GST the whole cluster can be used as an analogue ensemble. In contrast, for the variable PRES a low N a value is preferable.

5.3. Variation of the Analogue Size

The parameter k is used to define the number of records that constitute an analogue vector: that number is 2 k + 1 , with consecutive records in time separated by a Δ t time-step. Therefore, k also corresponds to half of the total number of time-steps within an analogue vector, and so k may also be referred to as the “half (time-)window size”.
All experiments discussed so far were conducted with k = 5 (analogue size = 11). In this section, the effects of varying k on the prediction errors are evaluated, while still considering the dependent stations scenario. Consistent with the results of the previous sections, the other relevant parameters were N c = 350 , N a c = 1 , N a = all cluster analogues.
Table 4 contains the results of varying the analogue size by varying k. The lowest errors (in bold) were obtained with a different set of k values for each variable. However, it is worth noting that k = 5 was shared by all variables and is therefore a good default value. This corresponds to a vector time span of 60 min between the first and the last sample. Therefore, one hour appears to be a reasonable time window to assess the correlation between the analogue candidates and the predictor and then make accurate predictions.
Another observation deserves an explanation: for k = 5 , the errors in Table 4 do not match the corresponding errors ( N a = ) in Table 3, although they were generated with the same parameters. This is because the values in these tables come from different runs of the K-means-based method, which is inherently non-deterministic, resulting in differences in the clusters formed in each run, which translates into differences in the errors produced.

5.4. Dependent vs. Independent Stations

Thus far, all results presented were obtained using the dependent station approach (Algorithm 4). In this section, a comparison is made with the results of Algorithm 3, which implements an independent-stations approach. All results refer to the same parameters: k = 5 , N c = 350 , N a c = 1 and N a = (all analogue candidates of the cluster used).
In Table 5, the best errors for each situation are shown in bold. For the independent approach, there is only one error of each type for each variable. For the dependent approach, two errors are given: the first is from row N a = of Table 3; the second (marked *) is from Table 4 and corresponds to the worst of the best values within the range of 10 2 given in Section 4.2. This is to show that even the worst of the best errors measured with the dependent approach are generally (with very few exceptions) smaller than the errors obtained with the independent approach.
From the results in Table 5—and in particular from the column of the combined error (CE)—it can be concluded that the use of the dependent method leads to a reduction of the prediction errors, even if this reduction was not very significant in some cases.

5.5. Prediction with Variables Different from the Predicted

In the previous sections, the data of the predictor stations corresponded to the data of the predicted variable in all predictions. For example, if the variable was ATMP, only the data from that variable were used to select the analogues. However, in a scenario where the stations used as predictors do not have historical records of the variable being predicted, the only alternative is to make predictions with a different variable, and the impact of this approach will vary depending on the correlation between the variables at stake.
Table 6 shows the prediction errors for each meteorological variable as a function of different predictor variables, keeping the same general parameters used so far: dependent approach, k = 5 , N c = 350 , N a c = 1 and N a = (all clusters analogue candidates used).
As can be seen in Table 6, the results for the variables PRES and ATMP are worse when a predictor variable other than the predicted one is used. However, for the WSPD and GST variables, the results obtained when one is used to predict the other are very similar to those for self-prediction, hinting that these two variables are highly correlated. These results suggest that the best analogues are obtained when predictors that are well correlated with the predicted variable are used.

6. Computational Performance

The classical AnEn method compares the predictor value with all historical values of the training period (step 1 of Figure 1) using the metrics of the Formulas (2) or (5). As discussed in Section 3.1, this requires a considerable amount of computation, especially if done purely sequentially. Furthermore, although it is possible to parallelize the analogues search [12], the longer the prediction period, the greater this effort will be.
A major advantage of determining analogues via K-means clustering is the potential for a dramatic reduction in overall computational time. As explained in Section 3.2, this reduction stems from the fact that clustering (which, since it is only done once, is still the most time-consuming phase in this approach), followed by the comparison of each cluster centroid with the predictor, is much faster than going through all possible analogues of the historical dataset and comparing them with the predictor.
This section describes the computational performance of implementations in R [22] of the classical and K-means-based AnEn methods that were specifically coded for this research. These implementations took advantage of built-in facilities in R for parallel execution, namely the parSapply function. Execution took place in a multicore environment provided by a Linux installation in a computer with two 16-core AMD EPYC 7351 CPUs (32 cores/64 threads in total) and 256 GB RAM.
For this performance study, the values chosen for the critical parameters of both AnEn methods were those previously shown to produce the lowest prediction errors. For the classical AnEN variant, these correspond to k = 2 and N a = 150 (recommended parameters in the paper [16]), while for the K-means variant, they correspond to those established in this paper (dependent approach, k = 5 , N c = 350 , N a c = 1 and N a = ).
The performance results are presented in Table 7 and Table 8. They unfold for the prediction of different variables, with different AnEn methods and a different number of CPU threads.
Table 7 allows a comparison of the best errors, total prediction times and the Speedup of the K-means variant against the classic method. The errors shown are limited to the Bias and SDE, which represent the systematic and random error, respectively. The results show that the K-means variant generally leads to more accurate predictions. Only for the variable PRES is this not the case, and even there the error differences are small. This is related to the discussion in Section 5.1 and Section 5.2, namely that the variable PRES requires a higher number of clusters due to its lack of variance for short time scales, so that each cluster characterizes a smaller sample, or alternatively, the number of subsets ( N a ) used as analogues in the cluster must be small.
With regard to the prediction times, these are similar for different variables when using the same AnEn variant. However, between the two AnEn variants, the prediction times are quite different and are up to two orders of magnitude apart, with the K-means variant being considerably faster (from 15- to 22-times faster) than the classical variant when using sequential processing ( n = 1 ), and keeping a lead when using parallel processing (although the K-means Speedup decreases with the number of CPU threads used).
In fact, the classical variant benefits greatly from the increase in the number of CPU threads used, since the computation of the similarity metric for each analogue candidate (step 1 in Figure 1) is independent of any other candidate and is thus performed in parallel. On the other hand, the K-means version does not benefit as much from the increased number of CPU threads, since only the similarity of centroids to the predictor can be computed in parallel, and their number is much smaller than the number of analogue candidates in the classical version.
Table 8 provides measures of the computational Speedup ( S n = T 1 / T n ) and Efficiency ( E n = S n / n ) as a function of the number (n) of CPU threads used. It can be observed that the classical approach has a higher Speedup and Efficiency for each n > 1 . Compared to the K-means approach, the classical approach has a Speedup that increases as n increases, and the loss of Efficiency is lower. As mentioned earlier, this method benefits greatly from the parallelization of analogue selection. However, to compete with the K-means variant in terms of prediction time, it needs a high number of CPU threads (as can be seen in Table 7, the K-means variant is still faster even with 64 CPU threads).
The results presented in this section show that the clustering-based variant of the AnEn method has better computational efficiency compared with the original version proposed by Monache [3,4]. Moreover, it is a reliable alternative that allows the reconstruction of missing data with errors smaller than or equal to the ones of the classical version.

7. Conclusions

The K-means AnEn variant is a worthy alternative to the classical AnEn method, with the same or better accuracy and much higher performance. This new variant facilitates using larger data sets and consequently solving larger problems in the various application domains of the AnEn method, such as hindcasting, forecasting and downscaling.
The numerical efficiency of this new AnEn variant depends on some of its most important parameters, namely the number of clusters to be formed ( N c ), the number of clusters effectively used ( N a c ), the number of analogues ( N a ) to be considered in each of the selected clusters and the time span (k) of each analog.
The way in which the data from multiple predictor stations are used to select the analogues (independently or with time dependence) and the use of the same or correlated predictor variables to make the predictions also have a measurable impact on their accuracy.
The experimental results showed that, for most of the meteorological variables considered, N c must be large enough for the subsets (vectors) included in each cluster to be sufficiently analogous to each other. Furthermore, it was shown that the use of a single best cluster ( N a c = 1 ) and all the subsets contained in it as analogues, together with an analogue time span of one hour ( k = 5 ), provided optimal or near-optimal accuracy.
As expected, the use of the same or highly correlated variables proved crucial to achieve the desired accuracy and the use of a dependent-stations approach proved beneficial (albeit slightly) for the same purpose.
Finally, from a computational-performance standpoint, the K-means-based AnEn approach has an undeniable advantage over the classical AnEn method: the prediction times of the former are much lower (up to two orders of magnitude) than those of the latter. However, the K-means variant doesn’t benefit from parallelization as much as the classical one. As the number of CPU threads used increases, the latter becomes more competitive: its Speedup scales better than that of the K-means variant, and the efficiency remains higher.
In the future, we want to investigate whether the application of machine-learning techniques brings advantages over the K-means-based variant proposed in this work.

Author Contributions

Conceptualization, C.B. and C.V.R.; methodology, C.B. and J.R.; software, J.R. and L.A.; validation, L.A., J.R. and C.V.R.; formal analysis, C.B. and C.V.R.; investigation, C.B. and L.A.; resources, J.R.; data curation, L.A. and C.V.R.; writing—original draft preparation, C.B.; writing—review and editing, J.R.; visualization, C.V.R.; supervision, C.B. and J.R.; project administration, C.B. and J.R.; funding acquisition, C.B. and J.R. All authors have read and agreed to the published version of the manuscript.

Funding

This work has been partially supported by FCT—Fundação para a Ciência e Tecnologia within the Project Scope: UIDB/05757/2020.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data used in this study are freely available from the United States National Data Buoy Center at https://www.ndbc.noaa.gov (accessed on 28 April 2022).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Lorenz, E.N. Atmospheric Predictability as Revealed by Naturally Occurring Analogues. J. Atmos. Sci. 1969, 26, 636–646. [Google Scholar] [CrossRef] [Green Version]
  2. Dool, H.M.V.D. A New Look at Weather Forecasting through Analogues. Mon. Weather Rev. 1989, 117, 2230–2247. [Google Scholar] [CrossRef] [Green Version]
  3. Monache, L.D.; Nipen, T.; Liu, Y.; Roux, G.; Stull, R. Kalman Filter and Analog Schemes to Postprocess Numerical Weather Predictions. Mon. Weather Rev. 2011, 139, 3554–3570. [Google Scholar] [CrossRef] [Green Version]
  4. Monache, L.D.; Eckel, F.A.; Rife, D.L.; Nagarajan, B.; Searight, K. Probabilistic Weather Prediction with an Analog Ensemble. Mon. Weather Rev. 2013, 141, 3498–3516. [Google Scholar] [CrossRef] [Green Version]
  5. Alessandrini, S.; Monache, L.D.; Sperati, S.; Nissen, J.N. A novel application of an analog ensemble for short-term wind power forecasting. Renew. Energy 2015, 76, 768–781. [Google Scholar] [CrossRef]
  6. Alessandrini, S.; Monache, L.D.; Sperati, S.; Cervone, G. An analog ensemble for short-term probabilistic solar power forecast. Appl. Energy 2015, 157, 95–110. [Google Scholar] [CrossRef] [Green Version]
  7. Cervone, G.; Clemente-Harding, L.; Alessandrini, S.; Delle Monache, L. Short-term photovoltaic power forecasting using Artificial Neural Networks and an Analog Ensemble. Renew. Energy 2017, 108, 274–286. [Google Scholar] [CrossRef] [Green Version]
  8. Rozoff, C.M.; Alessandrini, S. A Comparison between Analog Ensemble and Convolutional Neural Network Empirical-Statistical Downscaling Techniques for Reconstructing High-Resolution Near-Surface Wind. Energies 2022, 15, 1718. [Google Scholar] [CrossRef]
  9. Solomou, E.; Pappa, A.; Kioutsioukis, I.; Poupkou, A.; Liora, N.; Kontos, S.; Giannaros, C.; Melas, D. Analog ensemble technique to post-process WRF-CAMx ozone and particulate matter forecasts. Atmos. Environ. 2021, 256, 118439. [Google Scholar] [CrossRef]
  10. Yang, D.; van der Meer, D. Post-processing in solar forecasting: Ten overarching thinking tools. Renew. Sustain. Energy Rev. 2021, 140, 110735. [Google Scholar] [CrossRef]
  11. Vannitsem, S.; Bremnes, J.B.; Demaeyer, J.; Evans, G.R.; Flowerdew, J.; Hemri, S.; Lerch, S.; Roberts, N.; Theis, S.; Atencia, A.; et al. Statistical Postprocessing for Weather Forecasts: Review, Challenges, and Avenues in a Big Data World. Bull. Am. Meteorol. Soc. 2021, 102, E681–E699. [Google Scholar] [CrossRef]
  12. Hu, W.; Vento, D.; Su, S. Parallel Analog Ensemble—The Power of Weather Analogs. In Proceedings of the 2020 Improving Scientific Software Conference, Turin, Italy, 25–27 November 2020; pp. 1–14. [Google Scholar] [CrossRef]
  13. Hu, W.; Cervone, G.; Merzky, A.; Turilli, M.; Jha, S. A new hourly dataset for photovoltaic energy production for the continental USA. Data Brief 2022, 40, 107824. [Google Scholar] [CrossRef] [PubMed]
  14. Balsa, C.; Rodrigues, C.V.; Araújo, L.; Rufino, J. Hindcasting with Cluster-Based Analogues. In Communications in Computer and Information Science; Springer International Publishing: Berlin/Heidelberg, Germany, 2021; pp. 346–360. [Google Scholar] [CrossRef]
  15. Balsa, C.; Rodrigues, C.V.; Lopes, I.; Rufino, J. Using Analog Ensembles with Alternative Metrics for Hindcasting with Multistations. ParadigmPlus 2020, 1, 1–17. [Google Scholar] [CrossRef]
  16. Araújo, L.; Balsa, C.; Rodrigues, C.V.; Rufino, J. Parametric Study of the Analog Ensembles Algorithm with Clustering Methods for Hindcasting with Multistations. In Advances in Intelligent Systems and Computing; Springer International Publishing: Berlin/Heidelberg, Germany, 2021; pp. 544–559. [Google Scholar] [CrossRef]
  17. Chesneau, A.; Balsa, C.; Rodrigues, C.V.; Lopes, I.M. Hindcasting with multistations using analog ensembles. In Proceedings of the CEUR Workshop Proceedings (CEUR-WS), Copenhagen, Denmark, 30 March 2019; Volume 2486, pp. 215–229. [Google Scholar]
  18. Eldén, L. Matrix Methods in Data Mining and Pattern Recognition; SIAM: Philadelphia, PA, USA, 2007. [Google Scholar]
  19. Garey, M.R.; Johnson, D.S. Computers and Intractability—A Guide to the Theory of NP-Completeness; W. H. Freeman & Co.: New York, NY, USA, 1990. [Google Scholar]
  20. National Data Buoy Center. Available online: https://www.ndbc.noaa.gov/ (accessed on 15 April 2022).
  21. Chai, T.; Draxler, R.R. Root mean square error (RMSE) or mean absolute error (MAE)?—Arguments against avoiding RMSE in the literature. Geosci. Model Dev. 2014, 7, 1247–1250. [Google Scholar] [CrossRef] [Green Version]
  22. R Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2017. [Google Scholar]
Figure 1. Reconstruction of meteorological data with the analogues ensemble method ( k = 1 ).
Figure 1. Reconstruction of meteorological data with the analogues ensemble method ( k = 1 ).
Computation 10 00091 g001
Figure 2. Geolocation of the NDBC meteorological stations in Virginia (USA) [20].
Figure 2. Geolocation of the NDBC meteorological stations in Virginia (USA) [20].
Computation 10 00091 g002
Figure 3. Prediction errors as a function of the number of clusters ( N c ).
Figure 3. Prediction errors as a function of the number of clusters ( N c ).
Computation 10 00091 g003
Table 1. Meteorological stations.
Table 1. Meteorological stations.
StationCodeLocationRole
S 1 YKTV237°13 36 N 76°28 43 WPredicted
S 2 YKRV237°15 5 N 76°20 33 WPredictor
S 3 DOMV236°57 44 N 76°25 27 WPredictor
Table 2. Data availability in the meteorological stations.
Table 2. Data availability in the meteorological stations.
VariableMinMeanMax# NAAvailability
Station S 1 YKTV2
PRES [hPa]974.701017.341044.3020,15797.44%
ATMP [°C]−13.5016.0637.8026,10296.69%
WSPD [m/s]0.004.2623.8025,12796.81%
GST [m/s]0.005.4432.8025,15696.81%
Station S 2 YKRV2
PRES [hPa]972.601017.351043.9018,78597.62%
ATMP [°C]−12.8015.8636.3022,56697.14%
WSPD [m/s]0.005.9327.6024,79596.86%
GST [m/s]0.006.8839.6024,91096.84%
Station S 3 DOMV2
PRES [hPa]972.801017.771044.5019,99297.47%
ATMP [°C]−12.6016.1337.2025,19496.81%
WSPD [m/s]0.003.9124.3025,86696.72%
GST [m/s]0.005.2832.1025,89796.72%
Table 3. Prediction errors as a function of the number of analogues ( N a ).
Table 3. Prediction errors as a function of the number of analogues ( N a ).
N a BiasRMSEMAESDECEBiasRMSEMAESDECE
PRESATMP
500.5280.6790.5800.4272.214−0.0561.0490.7031.0472.855
1000.5620.7020.6130.4202.297−0.0551.0480.6961.0462.845
1500.5590.7030.6140.4272.303−0.0581.0470.6961.0452.846
2000.5610.7110.6170.4362.325−0.0491.0300.6911.0292.799
2500.5550.7060.6120.4362.309−0.0601.0480.6961.0462.850
3000.5570.7100.6140.4402.321−0.0521.0270.6911.0262.796
3500.5530.7050.6100.4362.304−0.0501.0300.6881.0282.796
4000.5650.7060.6180.4232.312−0.0531.0410.6971.0402.831
0.5610.7120.6180.4392.330−0.0371.0400.6951.0392.811
WSPDGST
50−0.1211.3771.0321.3723.902−0.3411.5791.1911.5424.653
100−0.0971.3661.0251.3623.850−0.3171.5641.1791.5324.592
150−0.0751.3621.0261.3603.823−0.2841.5711.1841.5454.584
200−0.0671.3571.0241.3553.803−0.2691.5631.1781.5394.549
250−0.0681.3651.0271.3643.824−0.2601.5541.1741.5324.520
300−0.0601.3521.0191.3513.782−0.2571.5611.1761.5404.534
350−0.0601.3581.0241.3573.799−0.2611.5571.1721.5354.525
400−0.0571.3641.0281.3623.811−0.2601.5531.1751.5314.519
−0.0591.3551.0241.3543.792−0.2631.5591.1721.5374.531
Table 4. Prediction errors as a function of the half time-window size (k).
Table 4. Prediction errors as a function of the half time-window size (k).
kBiasRMSEMAESDECEBiasRMSEMAESDECE
PRESATMP
10.5740.7170.6280.4302.349−0.0671.0430.6901.0412.841
20.5710.7120.6220.4252.330−0.0601.0380.6901.0362.824
30.5570.7040.6120.4302.303−0.0591.0460.6971.0442.846
40.5700.7130.6230.4282.334−0.0581.0490.6991.0472.853
50.550 *0.703 *0.607 *0.438 *2.298 *−0.053 *1.031 *0.685 *1.029 *2.798 *
60.5540.7080.6100.4412.313−0.0491.0390.6931.0382.819
70.5480.7010.6050.4372.291−0.0471.0480.6991.0472.841
80.5490.7050.6080.4422.304−0.0431.0450.7001.0442.832
90.5540.7080.6120.4412.315−0.0341.0440.7051.0432.826
100.5570.7130.6180.4452.333−0.0401.0420.6991.0412.822
WSPDGST
1−0.0631.3641.0281.3623.817−0.2581.5561.1721.5354.521
2−0.0601.3661.0291.3653.820−0.2581.558 *1.1691.537 *4.522
3−0.0601.3591.0261.3573.802 *−0.2571.5571.1721.5364.522
4−0.0621.3651.0291.3643.820−0.2561.558 *1.176 *1.537 *4.527 *
5−0.0551.361 *1.027 *1.359 *3.802 *−0.261 *1.5561.1711.5344.522
6−0.064 *1.3571.0251.3563.802 *−0.2581.5571.1711.5364.522
7−0.0611.3571.027 *1.3553.800−0.2661.5661.1801.5434.555
8−0.0661.3681.0321.3673.833−0.2581.5681.1861.5464.558
9−0.0621.3641.0301.3633.819−0.2511.5621.1801.5424.535
10−0.0711.3791.0381.3773.865−0.2571.5611.1761.5394.533
Table 5. Prediction errors with dependent and independent stations.
Table 5. Prediction errors with dependent and independent stations.
VariableDependencyBiasRMSEMAESDECE
PRESYes0.561, 0.550 *0.712, 0.703 *0.618, 0.607 *0.439, 0.438 *2.330, 2.298 *
No0.6530.7320.6850.3322.402
ATMPYes−0.037, −0.053 *1.040, 1.031 *0.695, 0.685 *1.039, 1.029 *2.811, 2.798 *
No−0.0511.0470.6851.0452.828
WSPDYes−0.059, −0.064 *1.355, 1361 *1.024, 1.027 *1.354, 1.359 *3.792, 3.811 *
No−0.0861.5581.1931.5554.392
GSTYes−0.263, −0.261 *1.559, 1.558 *1.172, 1.176 *1.537, 1.537 *4.531, 4.532 *
No−0.2781.6811.2901.6584.907
Table 6. Prediction errors for different variables used as predictors.
Table 6. Prediction errors for different variables used as predictors.
PredictedPredictorBiasRMSEMAESDECE
PRESPRES0.5640.7130.6190.4352.331
ATMP−0.6426.6515.2166.62019.129
WSPD−0.2106.5185.0856.51518.328
GST−0.1626.5905.1696.58818.509
ATMPPRES0.3437.9986.6217.99122.953
ATMP−0.0501.0380.6931.0372.818
WSPD−0.5588.7717.6708.75425.753
GST−0.5198.9167.8028.90126.138
WSPDPRES0.2222.5872.0672.5777.453
ATMP−0.0152.4431.9012.4436.802
WSPD−0.0591.3591.0241.3583.800
GST−0.0701.3271.0051.3253.727
GSTPRES0.0123.1962.5493.1968.953
ATMP−0.1983.0732.4023.0678.740
WSPD−0.2611.6091.2031.5884.661
GST−0.2571.5541.1711.5334.515
Table 7. Classical vs. K-means AnEn: prediction errors, computation times [s] and K-means Speedup.
Table 7. Classical vs. K-means AnEn: prediction errors, computation times [s] and K-means Speedup.
ErrorsNumber of CPU Threads (n) and Prediction Times ( T n )
VariableAnEnBiasSDE1248163264
PRESClassic0.2780.412206.0107.073.041.022.616.914.4
K-means0.5500.4389.06.05.04.03.73.43.7
Speedup22.8917.8314.610.256.114.973.89
ATMPClassic0.0001.070205.0111.052.039.021.217.015.1
K-means−0.0371.0299.06.05.04.23.43.33.7
Speedup22.7818.510.49.296.245.154.08
WSPDClassic−0.2062.064186.098.050.027.015.310.49.5
K-means−0.0591.35412.010.07.16.66.36.56.8
Speedup15.59.87.044.092.431.61.4
GSTClassic−0.5302.132200.098.053.027.015.310.98.3
K-means−0.2611.53713.011.08.06.66.76.16.7
Speedup15.388.916.634.092.281.791.24
Table 8. Classical AnEn vs. K-means AnEn: Speedup and Efficiency (%).
Table 8. Classical AnEn vs. K-means AnEn: Speedup and Efficiency (%).
Number of CPU Threads (n), Speedup ( S n ) and Efficiency ( E n )
VariableAnEnMeasure1248163264
PRESClassic S n 1.001.922.844.989.1212.2014.32
E n (%)100.096.271.062.357.038.122.4
K-means S n 1.001.471.912.152.382.592.38
E n (%)100.073.347.826.814.98.13.7
ATMPClassic S n 1.001.853.945.279.6712.0613.58
E n (%)100.092.798.465.960.537.721.2
K-means S n 1.001.511.932.122.622.702.41
E n (%)100.075.448.426.516.48.43.8
WSPDClassic S n 1.001.913.706.9312.1817.9219.62
E n (%)100.095.592.586.676.156.030.7
K-means S n 1.001.171.621.741.831.771.69
E n (%)100.058.740.521.811.45.52.6
GSTClassic S n 1.002.043.817.3113.0918.3824.13
E n (%)100.0102.095.291.481.857.437.7
K-means S n 1.001.201.591.911.882.071.88
E n (%)100.060.039.923.911.86.52.9
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Balsa, C.; Rodrigues, C.V.; Araújo, L.; Rufino, J. Cluster-Based Analogue Ensembles for Hindcasting with Multistations. Computation 2022, 10, 91. https://0-doi-org.brum.beds.ac.uk/10.3390/computation10060091

AMA Style

Balsa C, Rodrigues CV, Araújo L, Rufino J. Cluster-Based Analogue Ensembles for Hindcasting with Multistations. Computation. 2022; 10(6):91. https://0-doi-org.brum.beds.ac.uk/10.3390/computation10060091

Chicago/Turabian Style

Balsa, Carlos, Carlos Veiga Rodrigues, Leonardo Araújo, and José Rufino. 2022. "Cluster-Based Analogue Ensembles for Hindcasting with Multistations" Computation 10, no. 6: 91. https://0-doi-org.brum.beds.ac.uk/10.3390/computation10060091

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop