Scalable and Interpretable Forecasting of Hydrological Time-Series based on Variational Gaussian Processes

1. Introduction

Water resources represent one of the most fundamental natural resources. They play an indispensable role in sustaining life and human civilization, serving essential functions from facilitating agricultural irrigation to enabling power generation. Due to their importance, water resource forecasting models emerge for supporting operative activities such as flood mitigation, drought management, and hydropower optimization [40]. Besides, the water resource models offer a deeper understanding of hydrology, allowing ecological academic research [16,28] and informed decision-making [36]. Nowadays, streamflow forecasting emerges with force to mitigate climate change by supporting hydroelectric power generation [18], mainly in countries with geographically dispersed reservoir systems [1].

One country heavily dependent on hydroelectric power generation with a dispersed reservoir system is Colombia. The country projects a hydroelectric generation of more than 60% of the national demand by 2030 [8,41]. Colombia also conducts comprehensive studies on hydroelectric systems aiming at the established environmental sustainability objectives and policies for adopting renewable energies. Hence, studying these systems allows for a better understanding of their potential, efficient management, and contribution to environmental objectives, fostering the country’s energy autonomy. However, Colombia’s unique geographical and climatic characteristics introduce additional challenges in the operation and management of hydroelectric systems. The overall climatic diversity has caused high heterogeneity in hydrology, both temporally and spatially, presenting challenges in understanding and prediction [4]. Hence, hydrological models must explain the time-varying water resources and quantify their impact on hydroelectric planning and generation.

A first approach for addressing the prediction of water resources involves using deterministic physically-driven models, which establish the laws and natural interactions governing time series and thus present a high physical interpretation. However, these models enclose significant challenges as they require considerable effort to achieve acceptable results and depend on abundant information about specific parameters. Moreover, hydrological systems’ non-stationary and nonlinear characteristics lead to highly complex models with limited feasibility and intractable or complicated physical solutions [14,23,43]. Another source of complexity relies on the diverse time-series dynamics and data distributions among reservoirs that can be geographically distributed nationwide [27]. These intricacies necessitate the development of advanced, scalable forecasting models capable of adapting to the complex and stochastic nature of streamflow data [44].

Contrarily to physically-driven, statistical and machine learning models, mainly data-driven, allow for modeling inherent complexities in time series without involving physical processes [3]. Data-driven models have demonstrated higher performance than physically-driven in simulating streamflow within high-flow regimes, turning them into practical and accurate tools for hydrological regression [20].

Linear Autoregressive (LAR) models are the first choice among the data-driven models for time series forecasting due to their simplicity and interpretability. The linear autoregression predicts the forthcoming value of a target time series as a linear combination of its past values. Despite the acceptable performance on short-time hydrological forecasting, the augmented randomness and nonlinearity on longer prediction horizons overwhelm the LAR simplicity [35].

Long Short-Term Memory (LSTM) networks deal with LAR limitations by processing sequential data with nonlinear models of varying complexity. Unlike classical recurrent networks, LSTM addresses the gradient vanishing problem that often constrains the training from lengthy data sequences. The LSTM has been widely used in hydrological time series prediction, outperforming traditional artificial neural networks and conventional recurrent networks due to their ability to remember, forget, and update information [13,34]. For instance, an assembly of LSTM over other neural architectures enhanced individual deep models at hourly forecasting of the energy produced in the Songloulou hydropower plant.  [39]. However, the poor parameter interpretability turns these networks into black boxes [23]. Further, the hydrological power scheduling task involves strategic decision-making and action-planning processes that require quantifying the prediction uncertainty for alleviating operational and economic risks [12]. Unfortunately, the LSTM lacks uncertainty quantification, making it unable to identify outlying samples and warn of unreliable predictions [21].

To overcome the above LSTM constraints, the stochastic models provide predictive distributions instead of deterministic predictions [31]. In this regard, Gaussian Processes (GPs) are the stochastic model of choice for time series analysis. GPs associate a multivariate normal distribution to each time instant, appropriately representing the stochastic nature of the underlying phenomenon. Such models efficiently handle high-dimensional data projection through a kernel function that models temporal correlations and inherent nonlinearities [6]. An application in the energy-forecasting field was developed by [42] in which the authors used a Sparse Variational GP (SVGP) to describe day-ahead power wind and demonstrated its robustness. Another application for short-term electricity demand forecasting was proposed by [11] based on the SVGP and LSTM model. The results showed that the GP model achieves higher prediction accuracy than the LSTM model. Due to the extensive applicability of the Gaussian distributions in describing various phenomena and their ease of implementation, GPs have been successfully used in hydrological time series prediction, yielding satisfactory results [29,37,38].

Furthermore, real-world scenarios demand the simultaneous prediction of multiple time series, an approach known as multi-output forecasting. A straightforward multi-output model extends the single-output model to generate various predictions simultaneously instead of training several output-wise models. Since all target variables influence the parameter training of such a multi-output approach, the outputs borrow statistical dependencies. For instance, extending nonlinear regression models for incorporating layers dedicated to each task improved the simultaneous forecasting of multiple time series [30]. Another multi-output model implemented an auxiliary task that supported the prediction of the main output while providing an additional regularization [24]. Also, the integrated forecasting approach of wind, photovoltaic, and hydropower energy improves modeling efficiency compared to independent task forecasting [10]. In the same spirit, the GP model extension to multi-output scenarios allowed approximating the correlation among target outputs through kernel functions [45].

Nonetheless, the MOGP approach increases its tunable and trainable parameters, intensifying the complexity of modeling target correlations [2,25]. Particularly for hydrological applications, such complex models generally outperform simpler ones at the cost of reduced stability due to an over-parameterized structure for extracting the time series information. Then, the over-parameterization poses challenges for reliable parameter estimation and validation in multi-output forecasting of hydrological time series [9].

Considering the previous challenges, this work develops a methodology for forecasting multiple hydrological reservoir streamflows based on Sparse Variational Gaussian Processes (SVGPs). The proposed SVGP-based forecasting includes a kernel-based covariance matrix to exploit target dependencies and nonlinearities. Further, the variational treatment yields a set of inducing points shared by all outputs that reduce the model complexity compared to conventional multi-output Gaussian Processes. The kernel-based covariance and shared inducing points endow the SVGP with interpretability thanks to the automatic relevance determination criterion and latent data distribution modeling. This work also contrasts the SVGP against baseline Linear Autoregressive and Long-Short Term Memory models in the simultaneous forecasting of 23 streamflow time series from Colombian reservoirs daily recorded over twelve years. The three contrasted approaches operate within a multi-output forecast framework from one to thirty-day horizons. The attained results prove the advantage of the proposed SVGP-based approach over LAR and LSTM by reducing the mean-squared prediction error in a two-year testing period with three additional benefits: The provided predictive distribution to handle forecast uncertainty, the inherent modeling of reservoir relationships through kernel functions and shared inducing points, and the reduction of model complexity through the variational strategy and independent kernel.

This work considers the following agenda for developing the proposed methodology: Methodology details the theoretical foundations for the SVGP-based forecasting methodology. Results develops the hyperparameter tuning, the interpretability analysis, and the performance assessment. Conclusions concludes the work with the findings and future research directions.

2. Mathematical Framework

2.1. Sparse Gaussian Process Regression Framework

The Gaussian Process (GP) model offers an elegant framework for formulating non-parametric and probabilistic fashion regression tasks, yielding a point estimate and capturing the associated uncertainty through a Gaussian distribution. This model assigns a multivariate Gaussian distribution to every finite set of output points, thereby capitalizing on inherent correlations. By taking advantage of the Gaussian distribution, it becomes feasible to establish a prior distribution over the target function, which can then be conditioned on an observation dataset to derive the posterior distribution, effectively integrating the prior distribution knowledge and the data information.

Consider a dataset denoted as $\mathcal{D}={\{{\mathit{x}}_n,{\mathit{y}}_n\}}_{n=1}^N$ comprising N input-output pairs, where ${\mathit{x}}_n\in \mathit{\chi }$ and $\mathit{\chi }\subset {\mathbb{R}}^L$ represents the L-dimensional input space. The target function is a vector-valued function where ${\mathit{y}}_n\in {\mathbb{R}}^D$ comprises the observations of all outputs (tasks) at the same input step. In the context of the GP framework, the dataset $\mathcal{D}$ supports the learning of a random mapping function $\mathit{f}(·)$ that captures the relationship between ${\mathit{x}}_n$ and ${\mathit{y}}_n$, and follows a zero mean prior distribution:

$$\mathit{f}\left(\mathit{x}\right)\sim \mathcal{GP}\left(\mathit{0},\mathit{k}(\mathit{x},{\mathit{x}}^{\prime }\mid \mathit{\theta })\right)$$

(1)

where $\mathit{f}:\mathit{\chi }\to {\mathbb{R}}^D$ represents the vector-valued function responsible for mapping the input space, and $\mathit{k}:\mathit{\chi }\times \mathit{\chi }\to {\mathbb{R}}^{D\times D}$ symbolizes the across-tasks covariance matrix function, parameterized by vector $\mathit{\theta }$. Adding i.i.d zero mean Gaussian noise $\mathit{ϵ}$ with covariance matrix ${\sigma }_{ϵ}^2{\mathit{I}}_D$ as a regularization term, the model prediction given by mogpnotation turns into ${\mathit{y}}_n=\mathit{f}\left({\mathit{x}}_n\right)+\mathit{ϵ}$, corresponding to a noise observation of the latent variable $\mathit{f}$.

The matrix function, $\mathit{k}$, can be systematically described in terms of its components. To this end, consider the expansion of the input space as ${\mathit{\chi }}_{\mathrm{ext}}=\mathit{\chi }\times \{1,\cdots ,D\}$. The above allows for a reconfiguration of the dataset, transforming it into an assembly of $ND$ input-output pairs denoted as $\mathcal{D}={\{{{\mathit{x}}_{\mathrm{ext}}}_n,y_n\}}_{n=1}^{ND}=\{\mathit{X},\mathit{y}\}$, where ${{\mathit{x}}_{\mathrm{ext}}}_n\in {\mathit{\chi }}_{\mathrm{ext}}$ and $y_n$ stands for a single task observation.

With this extended framework, the vector-valued prior presented in mogpnotation can be perceived as a single output model, spanning over an extended input space, as follows [2]:

$$f\left({\mathit{x}}_{\mathrm{ext}}\right)\sim \mathcal{GP}\left(0,k({\mathit{x}}_{\mathrm{ext}},{\mathit{x}}_{\mathrm{ext}}^{\prime }\mid \mathit{\theta })\right)$$

(2)

being $k:{\mathit{\chi }}_{\mathrm{ext}}\times {\mathit{\chi }}_{\mathrm{ext}}\to \mathbb{R}$ the scalar kernel function. Let be ${\mathit{X}}_{*}={\left\{{{\mathit{x}}_{\mathrm{ext}}}_n\right\}}_{n=1}^{N_{*}D}$ a collection of test point with test output vector ${\mathit{f}}_{*}\in {\mathbb{R}}^{N_{*}D}$ as the same way as $\mathit{y}$. The joint Gaussian distribution over the above vector is specified as follows:

$$\begin{array}{c}\hfill \left[\begin{array}{c}\mathit{y}\\ {\mathit{f}}_{*}\end{array}\right]\sim \mathcal{N}\left(\mathit{0},\left[\begin{array}{cc}{\mathit{K}}_y& {\mathit{K}}_{*}\\ {\mathit{K}}_{*}^{\top }& {\mathit{K}}_{**}\end{array}\right]\right)\end{array}$$

(3)

Here, the vector ${\mathit{f}}_{*}$ is associated with a prior distribution, and ${\mathit{K}}_y=\mathit{K}+{\sigma }_{ϵ}^2{\mathit{I}}_{ND}$ where $\mathit{K}$ is formed by evaluation of the scalar covariance function at all pairs of input vectors in $\mathit{X}$, and ${\mathit{K}}_{**}$ and ${\mathit{K}}_{*}$ are formed evaluating the covariance function across all pairs of test inputs in ${\mathit{X}}_{*}$ and test-train inputs respectively. This notation allows deriving the conditional distribution named posterior as follows:

$${\mathit{f}}_{*}|{\mathit{X}}_{*},\mathcal{D}\sim \mathcal{N}(\overline{{\mathit{f}}_{*}},\mathrm{cov}\left({\mathit{f}}_{*}\right))$$

(4)

with

$$\overline{{\mathit{f}}_{*}}={\mathit{K}}_{*}^{\top }{\mathit{K}}_y^{-1}\mathit{y}$$

(5)

$$\mathrm{cov}\left({\mathit{f}}_{*}\right)={\mathit{K}}_{**}-{\mathit{K}}_{*}^{\top }{\mathit{K}}_y^{-1}{\mathit{K}}_{*}$$

(6)

Notice from sogpposteriormean,sogpposteriorcov that the mean function lies near the observations $\mathit{y}$ and the variance decreases where test points are close to train points. In contrast, the posterior distribution tends to the prior for test instances far from the training ones, merging the data and prior knowledge. Furthermore, the prediction performance archived by the conditional distribution depends on the selected parameter set $\mathit{\theta }$ and observation noise variance ${\sigma }_{ϵ}^2$. Both parameters are calculated by maximizing the marginal log-likelihood from sogppriormatrix where $p\left(\mathit{y}\right)=\mathcal{N}(\mathit{y}\mid \mathit{0},{\mathit{K}}_y)$ as follows:

$$\begin{array}{cc}\hfill \{{\mathit{\theta }}_{\mathrm{opt}},{{\sigma }_{ϵ}^2}_{\mathrm{opt}}\}& =\underset{\mathit{\theta },{\sigma }_{ϵ}^2}{argmax}lnp\left(\mathit{y}\right)\hfill \\ \hfill & =\underset{\mathit{\theta },{\sigma }_{ϵ}^2}{argmin}\frac{1}{2}{\mathit{y}}^{\top }{\mathit{K}}_y^{-1}\mathit{y}+\frac{1}{2}ln\mid {\mathit{K}}_y\mid +\frac{ND}{2}ln2\pi \hfill \end{array}$$

(7)

The optimization problem in sogpnlmlopt is conventionally addressed using a gradient-based optimizer [33]. However, the problem holds cubic complexity of $\mathcal{O}\left(N^3{D}^3\right)$ and storage demand of $\mathcal{O}\left(N^2{D}^2\right)$ due to inverting the matrix ${\mathit{K}}_y$, a procedure usually accomplished through Cholesky decomposition. That complexity becomes impractical for datasets containing a few thousand samples [15].

2.2. Sparse Variational Gaussian Process

To alleviate the complexity, the sparse variational strategy defines a set of trainable points $\mathit{Z}={\{{\mathit{z}}_m\in {\mathit{\chi }}_{\mathrm{ext}}\}}_{m=1}^{MD}$, termed inducing points, such that $M\ll N$, and laying in the extended input space ${\mathit{\chi }}_{\mathrm{ext}}$. Such a treatment also determines a respective inducing variable vector $\mathit{u}={\left\{f\left({\mathit{z}}_m\right)\right\}}_{m=1}^{MD}\in {\mathbb{R}}^{MD}$, playing a similar role to the training output vector $\mathit{f}={\left\{f\left({\mathit{x}}_n\right)\right\}}_{n=1}^{ND}$ that holds computationally expensive marginal distribution $p\left(\mathit{f}\right)={\mathbb{E}}_{p\left(\mathit{u}\right)}\left\{p(\mathit{f}\mid \mathit{u})\right\}$. The variational approach approximates the posterior $p\left(\mathit{f}\right|\mathit{y})$ through $q\left(\mathit{f}\right)$ by averaging the conditional distribution of $\mathit{f}$ w.r.t. the inducing variable vector model $q\left(\mathit{u}\right)$, yielding a Sparse Variational Gaussian Process:

$$q\left(\mathit{f}\right):={\mathbb{E}}_{q\left(\mathit{u}\right)}\left\{p(\mathit{f}\mid \mathit{u})\right\}$$

(8)

Since the prior output distribution $p\left(\mathit{f}\right)$ is Gaussian, $q\left(\mathit{u}\right)$ becomes also Gaussian with mean vector $\mathit{m}\in {\mathbb{R}}^{MD}$ and covariance matrix $\mathit{S}\in {\mathbb{R}}^{MD\times MD}$, so that the approximating distribution $q\left(\mathit{f}\right)$ has a closed form in prioraproximation given by:

$$q\left(\mathit{f}\right)=\mathcal{N}\left(\mathit{f}\mid \mathit{A}\mathit{m},\mathit{K}+\mathit{A}(\mathit{S}-{\mathit{K}}_{MM}){\mathit{A}}^{\top }\right)$$

(9)

with $\mathit{A}={\mathit{K}}_M{\mathit{K}}_{MM}^{-1}$, ${\mathit{K}}_M\in {\mathbb{R}}^{ND\times MD}$ as the scalar kernel function in sogpnotation evaluated for every combination of inducing and training points. ${\mathit{K}}_{MM}\in {\mathbb{R}}^{MD\times MD}$ corresponds to the kernel function evaluated at all inducing point pairs. As a result, the tuning for parameters of $q\left(\mathit{u}\right)$, $\mathit{\theta }$, and the observation noise variance ${\sigma }_{ϵ}^2$ becomes a stochastic framework, optimizing the manageable lower bound of marginal log-likelihood within the cost function in sogpnlmlopt as follows:

$$lnp\left(\mathit{y}\right)\ge \sum _{n=1}^{ND}{\mathbb{E}}_{q\left(f_n\right)}\left[lnp(y_n\mid f\left({\mathit{x}}_n\right))\right]-\mathrm{KL}\left[q\left(\mathit{u}\right)\Vert p\left(\mathit{u}\right)\right]$$

(10)

being KL the Kullback-Leibler divergence measuring the relative entropy between the distributions $q\left(\mathit{u}\right)$ and $p\left(\mathit{u}\right)$ and also playing the role of regularization term. This approach provides a significant advantage by reducing the complexity of the GP model to $\mathcal{O}\left(N{M}^2{D}^3\right)$ through the inversion of a smaller matrix ${\mathit{K}}_{MM}$. Consequently, the model efficiently manages larger datasets at a lower overall cost.

2.3. Model Setup

The proposed approach factorizes the GP scalar covariance function in sogpnotation on two kernel functions: $k_{\mathit{\chi }}$ models input correlations through relative distances [45], and $k_D$ models tasks pair-wise correlation as following in separablekernel:

$$\begin{array}{cc}\hfill k\left({\mathit{x}}_{\mathrm{ext}},{\mathit{x}}_{\mathrm{ext}}^{\prime }\mid \mathit{\theta }\right)& =k\left((\mathit{x},d),({\mathit{x}}^{\prime },d^{\prime })\mid \mathit{\theta }\right)\hfill \end{array}$$

$$\begin{array}{cc}\hfill & =k_{\mathit{\chi }}\left(\mathit{x},{\mathit{x}}^{\prime }\mid {\Theta }_d\right)k_D\left(d,d^{\prime }\mid {\sigma }_d\right)\hfill \end{array}$$

(11a)

$$\begin{array}{cc}\hfill k_{\mathit{\chi }}\left(\mathit{x},{\mathit{x}}^{\prime }\mid {\Theta }_d\right)& =\exp \left(-\frac{1}{2}{(\mathit{x}-{\mathit{x}}^{\prime })}^{\top }{\Theta }_d^{-2}(\mathit{x}-{\mathit{x}}^{\prime })\right)\hfill \end{array}$$

(11b)

$$\begin{array}{cc}\hfill k_D\left(d,d^{\prime }\mid {\sigma }_d^2\right)& ={\sigma }_d^2{\delta }_{d,d^{\prime }}\hfill \end{array}$$

(11c)

where ${\delta }_{d,d^{\prime }}$ is the Kronecker delta function between task indices, ${\sigma }_d^2$ is the output scale corresponding with task d and the diagonal matrix ${\Theta }_d=\mathrm{diag}{\left\{{\Delta }_{dl}\right\}}_{l=1}^L$ gathers the positive lengthscale factors ${\Delta }_{dl}\in {\mathbb{R}}^{++}$ from each input dimension corresponding with task d and thus, the trainable covariance parameters becomes into $\mathit{\theta }={\{{\sigma }_d^2,{\Theta }_d\}}_{d=1}^D$. Notice that the input kernel $k_{\mathit{\chi }}$ implements the widely used squared-exponential kernel function in squaredexponentialkernel, allowing a smooth data mapping. In contrast, the task kernel $k_D$ given in kroneckerkernel does not model the correlations between tasks to avoid making the model unstable due to a high amount of hyperparameter to train and making the matrix ${\mathit{K}}_{MM}$ diagonal by blocks, reducing the computational complexity even more to $\mathcal{O}\left(N{M}^{2}D\right)$. Although the correlations are not modeled by the output kernel $k_D$, because the autoregressive model uses as inputs the values of the time series in previous steps, the lengthscale values of the kernel $k_{\mathit{\chi }}$ and shared inducing points set may support the exploration of dependencies between tasks with fewer parameters.

3. Results and Discussions

3.1. Dataset Collection

The hydrological forecasting task for validating the proposed Sparse Variational Gaussian Process (SVGP) regressor considers time series data of streamflow contributions from 23 Colombian reservoirs. These streamflow contributions were recorded daily from January 1st, 2010, to February 28th, 2022, resulting in a dataset of 4442 samples. It is worth noting that while these contributions represent volumetric values, they are reported in kilowatt-hours (kWh) by the hydroelectric power plants, which is a more practical unit for their daily operational considerations. Figure 1 locates the reservoirs on a hydrological map for contextualizing the task. Note that reservoirs distribute along the mountains to profit from the streamflow produced by every precipitation.

Figure 2 depicts the reservoir-wise data distribution. Firstly, note the wide span of streamflow contributions among reservoirs due to differing power generation capacities. Secondly, most reservoirs hold a non-Gaussian distribution, including platykurtic (such as reservoir C), negatively skewed (reservoir A), positively skewed (reservoir B), long-tailed (reservoir F), and bimodal (reservoir O). The factors leading to the above varying distributions include the atypical rainy days, prolonged dry seasons prevalent in Colombian weather conditions, and operational decisions for enforcing water conservation or maximizing the generation. Further, reservoir Q, which consistently reported zero streamflow, allows for assessing the forecasting performance in trivial cases. For validation purposes, the last two years of data assess the SVGP performance as the testing subset, while the first ten years become the training subset with $N=3554$ samples.

3.2. Hyperparameter Tuning and Model Interpretability

Given the Colombian streamflow dataset, the hydrological forecasting task corresponds to predicting the $D=23$ reservoir contributions at time instant $n+H$ from the $L=23$ contributions at day n. H stands for the forecasting horizon, including short (up to a week ahead $H\in \{1,2,3,4,5,6,7\}$) and medium-term (up to one month ahead $H\in \{9,14,21,30\}$ ). To this end, the SVGP demands a fixed number of inducing points, M. Figure 3 presents the grid search results for tuning M by minimizing the MSE in a ten-fold cross-validation on the training subset at the fixed one-day-ahead forecasting horizon ($H=1$). Results evidence that the median MSE generally decreases when increasing the number of inducing points, up to $M=16$, after which the error increases. The standard deviation also decreases within the range $M\in \{1,\cdots ,16\}$, suggesting a consistency improvement. The above is because a few inducing points are more sensitive to initialization and fail to capture the intrinsic phenomenon distribution, thereby underfitting the regression task. On the contrary, too many inducing points monotonically increase the median and standard deviation on the error, indicating a model overfitting. Consequently, the hyperparameter tuning designates $M=16$ inducing points as the most accurate and consistent within the grid search.

Using the tuned number of inducing points, Figure 4 presents the trained reservoir-wise output scales ${\sigma }_d^2$ and the estimated noise variance ${\sigma }_{ϵ}^2$ at each horizon h. Overall, the longer the horizon, the smaller the output scale for all reservoirs. Such a fact suggests a reduction of time correlation for long horizons, making the SVGP approximate the prior distribution with a time-uncorrelated function. Nonetheless, for reservoirs P and V, an output scale increment at long horizons follows the initial reduction, implying the existence of relevant periodicity components within such time series. Lastly, reservoir Q consistently renders a zero scale along horizons for enforcing the expected zero streamflow.

Regarding the noise variance, the heightened unpredictability in the output data at long horizons yields an increment in ${\sigma }_{ϵ}^2$, which also supports the positive definiteness of the covariance matrix ${\mathit{K}}_y$ in sogpposteriormean. For such horizons, the model uncertainty mainly depends on the noise variance rather than the output scale. Therefore, the SVGP model exploits the data correlations at short horizons, yielding a low uncertainty while understanding the lack of information for long-horizon forecasting as a zero-mean isotropic Gaussian posterior.

The lengthscale analysis relies on its capacity to deactivate nonessential features on a kernel machine according to the Automatic Relevance Determination (ARD) criterion [7,26]. Hence, a small lengthscale value implies a high relevance for predicting the future streamflow of a reservoir from the past of another [5]. Figure 5 depicts the trained lengthscales for squaredexponentialkernel for each input feature (columns) to each output task (rows) for prediction horizons of one day (left), two weeks (center), and one month (right). At first glance, reservoir Q exhibits notably low lengthscales, either for predicting other reservoirs (Q-th column) or as the target reservoir (Q-th row), misinterpreted as highly relevant according to ARD. As the target reservoir, the zero streamflow during the observation period tunes to zero the output scale, turning off the corresponding rows in ${\mathit{K}}_{*}$ (in sogpposteriormean) and the influence of every reservoir on the posterior mean for the reservoir Q despite its length scale. As the predicting reservoir, reservoir Q holds N zero terms within the vector $\mathit{y}$ in sogpposteriormean, yielding a zero contribution to build the mean posterior of every other reservoir. In either case, the SVGP identifies the influence lack of the reservoir Q and fixes the corresponding lengthscales (rows and columns) at the default value of ${\Delta }_{dl}=1$.

The shortest horizon (one day ahead prediction, $H=1$) distributes the smallest lengthscale values over the main diagonal contrary to the high off-diagonal values. This result proves the relevance of short-term memory components within each streamflow. At a longer horizon of fourteen days ($H=14$), the lengthscales over the main diagonal become higher, losing relevance, and several off-diagonal lower, gaining relevance, than at a one-day horizon. Such a change implies that a reservoir provides less information to predict itself at longer horizons, leading the model to look for knowledge on other reservoirs. The longest horizon ($H=30$) illustrates even less relevant lengthscales, indicating that the SVGP restrains the number of relevant features. Thus, the lengthscale variations across different prediction horizons reveal the dynamic interplay between reservoirs and highlight the Gaussian Process’s ability to select relevant features adaptively.

For visual interpretation of SVGP inner workings, Figure 6 maps in 2D the reservoir data along with the optimized inducing points using the t-distributed Stochastic Neighbor Embedding (t-SNE) technique for each target reservoir. The input-space similarity measure demanded by t-SNE was tailored to the kernel in squaredexponentialkernel with the tuned lengthscales for each task, guaranteeing the mapping of the same latent space designed by the SVGP. The filled contours outline the density distribution of training data, while the color-numerated scatter locates the optimized inducing points. Since the forecasting tasks share the inducing points, not the kernel parameters, each target reservoir holds a different t-SNE mapping.

The 2D plots exhibit multi-modal distributions for all tasks, indicating clustered information and varying dynamics within a reservoir. Hence, the SVGP trades off the location of the inducing points between the within-task centroids and the among-task relationships. For instance, the seventh inducing point emerges near a mode of the four distribution plots, representing global information shared by most tasks. The locations of points 12 and 16 also pose a remarkable insight. On the one hand, the former mainly represents one of the modes of reservoirs E (Figure 6a) and U (Figure 6d) simultaneously, despite being located far from any cluster of I and L. On the other hand, the MOVGP devotes the latter point for the most prominent mode of reservoirs I (Figure 6b) and L (Figure 6c), but missing E and U. Nonetheless, those two inducing points (12-th and 16-th) provide complementary data representativity for improving the performance of a few related tasks. Therefore, the shared inducing points allow capturing the task-wise, group-wise, and global information about the streamflow dynamics.

3.3. Performance Analysis

For performance analysis, the proposed SVGP-based forecasting is contrasted against the straightforward Linear AutoRegressive model (LAR) and the nonlinear Long-Short-Term Memory network (LSTM) [17,19,22,32,34]. Figure 7 illustrates the testing time series of streamflow contributions and their corresponding one-day ahead predictions generated by the contrasted models for three reservoirs with different forecasting complexities. The red dots correspond to test data, and the yellow and green lines denote the LAR and LSTM predictions respectively. The blue line and region represent the mean and 95% confidence interval of the SVGP predictive distribution.

Reservoir U in Figure 7a presents the challenge of changing dynamics within the time series, from a steady low variance to a rapidly evolving streamflow. For such a reservoir, SVGP outperforms by better adapting the forecasted dynamics. Reservoir E in Figure 7b portrays a scenario with reduced smoothness with a few abruptly changing segments. Even for those segments, the proposed model tracks the actual curve more closely than LAR and LSTM. Further, the SVGP effectively explains the existence of sudden shifts through its intrinsic confidence intervals. In contrast, reservoir I in Figure 7c exhibits a heavily varying streamflow. Despite the forecasting for the three models diverging from the actual peak values within the time series, the predictive distribution provided by the SVGP interprets them as outliers. In general, the SVGP outstands the most adaptative and suitable model to exploit the changing within- and between-reservoir time series, explaining their stochastic nature through the inherent predictive distribution.

Finally, Figure 8 displays the overall Mean Squared Error (MSE) scores achieved by the baseline approaches of LAR and LSTM and the proposed SVGP-based streamflow forecasting over a testing period of two years. Each row represents one target reservoir, while the columns indicate the forecasting horizon. The lower plot indicates the grand-average MSE along the horizon. The colorbar on the right maps the resulting MSE to a logarithmic scale to favor the visual interpretation.

As expected, the longer the horizon, the larger the error for all approaches, implying more challenging future forecasting. However, the prediction difficulty varies among reservoirs, with D, H, I, and N as the most demanding, suggesting a lack of information from external factors. In contrast, reservoirs A, K, and U exhibit a less complex streamflow dynamic due to a lower MSE. A remarkable observation pertains to reservoir Q, for which the SVGP model sustains a zero error due to the null contribution of its corresponding time series in contrast to other models, which obtain a residual error.

From one to seven days horizons, the proposed SVGP outperforms LAR and LSTM in reservoirs A, C, E, F, and K. Hence, the SVGP effectively captures short-term dependencies in the data. Particularly on the reservoir E, the SVGP considerably reduces the error of LSTM, indicating a benefit in capturing short-term dynamics compared to a nonlinear model. As the prediction horizon extends to 14 to 30 days, the SVGP performance remains competitive since its MSE increases slower than LAR and LSTM along the horizon, especially in reservoirs like L, T, and V. Such a fact can also be seen in the grand average MSE where SVGP error remains considerably slower than the baseline models. Concludingly, the Gaussian Process excels in all tested horizons mainly due to the capabilities of nonlinear mapping the data into an infinitely high-dimensional feature space and task-wise adaptation to varying dynamics.

4. Concluding Remarks

This work develops a scalable and interpretable streamflow forecasting methodology based on Sparse Variational Gaussian Processes (SVGP). Compared to the conventional GPs, the proposed methodology diminishes the computational complexity from cubic to linear regarding the number of samples, becoming more scalable to forecasting tasks on large datasets. In this sense, the proposed SVGP focuses on predicting 23 Colombia’s reservoir streamflow contributions, geographically distributed nationwide. Due to the widespread locations, time-series dynamics and data distribution vary among reservoirs.

Regarding the hyperparameters, the optimal number of inducing points performs as an inherent regularization by encoding the most information from the data while avoiding overfitting. Such an optimal hyperparameter allows interpretation of the Gaussian Process model in terms of the inducing points and the kernel parameters. On the one hand, the t-SNE-based distribution plots reveal that the proposed model strategically places the shared inducing points to capture task-wise, group-specific, and global dynamics despite the complex, multi-modal nature of reservoir data, trading off between capturing the reservoir-wise unique characteristics and the between-reservoir dependencies, thereby improving the overall performance in streamflow forecasting. On the other hand, the trained lengthscales prove that the GP model successfully adjusts its focus to prediction horizon changes between short-term memory components (using the past data of the reservoir itself) and long-term relationships (exploiting interactions between different reservoirs). This adaptability makes the Gaussian Process model robust and versatile for multiple-output forecasting tasks with varying time features.

The performance analysis evidences the advantage of the SVGP-based model over the baseline Linear Autoregressive (LAR) model and Long Short-Term Memory network (LSTM) for streamflow prediction in three main aspects. The performance analysis evidences the advantage of the proposed SVGP model over the baseline Linear Autoregressive (LAR) model and nonlinear Long Short-Term Memory network (LSTM) for streamflow prediction in three main aspects. Firstly, the proposed approach suitably adapts to changing dynamics within and between reservoirs, even on rapidly changing time series or abrupt shifts. Secondly, the Bayesian scheme makes the proposed methodology predict the streamflow as the mean of an intrinsic posterior distribution and provides informative confidence intervals through its posterior variance. Thirdly, the SVGP better copes with the enhanced forecasting complexity for long horizons than baseline approaches, as proved by the slower error growth. Overall, the proposed SVGP-based methodology efficiently captures and adapts to the complex and stochastic nature of the streamflow contribution time series, enhancing the forecasting task on multiple horizons.

For pursuing the research on SVGP-based hydrological forecasting, three directions are advised: The first research direction aims to integrate SVGP-based streamflow forecasting into the scheduling algorithms of thermoelectric power plants. Such an integrated approach will support a resilient and efficient energy system by optimizing the operation of hydroelectric and thermoelectric resources. Another future work is related to the fact that the reservoirs become operational at different dates. Then, the SVGP-based forecasting must focus on maintaining a suitable performance on missing data scenarios. One last research in the artificial intelligence field explores how to provide non-Gaussian likelihoods following better fitting the constraints of the output distribution, e.g., non-negative streamflow data. Such an approach can be extended to a chained GP framework, allowing each likelihood parameter to follow a GP and enhancing the model flexibility to handle the complex dynamic of hydrological time series.