AI-based techniques for multi-step streamflow forecasts: application for multi-objective reservoir operation optimization and performance assessment

Streamflow forecasts are traditionally effective in mitigating water scarcity and flood defense. This study developed an artificial intelligence (AI)-based management methodology that integrated multi-step streamflow forecasts and multi-objective reservoir operation optimization for water resource allocation. Following the methodology, we aimed to assess forecast quality and forecast-informed reservoir operation performance together due to the influence of inflow forecast uncertainty. Varying combinations of climate and hydrological variables were input into three AI-based models, namely a long short-term memory (LSTM), a gated recurrent unit (GRU), and a least-squares support vector machine (LSSVM), to forecast short-term streamflow. Based on three deterministic forecasts, the stochastic inflow scenarios were further developed using Bayesian model averaging (BMA) for quantifying uncertainty. The forecasting scheme was further coupled with a multi-reservoir optimization model, and the multi-objective programming was solved using the parameterized multi-objective robust decision-making (MORDM) approach. The AI-based management framework was applied and demonstrated over a multi-reservoir system (25 reservoirs) in the Zhoushan Islands, China. Three main conclusions were drawn from this study: (1) GRU and LSTM performed equally well on streamflow forecasts, and GRU might be the preferred method over LSTM, given that it had simpler structures and less modeling time; (2) higher forecast performance could lead to improved reservoir operation, while uncertain forecasts were more valuable than deterministic forecasts, regarding two performance metrics, i.e., water supply reliability and operating costs; (3) the relationship between the forecast horizon and reservoir operation was complex and depended on the operating configurations (forecast quality and uncertainty) and performance measures. This study reinforces the potential of an AI-based stochastic streamflow forecasting scheme to seek robust strategies under uncertainty.

1 Introduction

Multi-step streamflow forecast is of great importance for reservoir operations to determine optimal water allocations considering the current use and the carry-over storage for mitigating water scarcity risk in the future (Guo et al., 2018; Zhao et al., 2019). Previous studies have identified that real-time reservoir operations are influenced by multiple uncertainties (Xu et al., 2020), among which inflow forecast uncertainty has been determined as the primary source, resulting in the risk of water shortage when the forecast inflow overestimates the actual inflow. Ensemble forecasting techniques are commonly used to characterize various uncertainties in streamflow forecasts. According to comparative analysis for various probabilistic forecasting techniques (Nott et al., 2012; Fang et al., 2018a; Zhai and Chen, 2018; Zhou et al., 2020b), Bayesian model averaging (BMA) (Hoeting et al., 1999) has been found to be an effective and most commonly used method to evaluate uncertainty and thus can be used in streamflow forecast.

Any ensemble forecast approach relies upon model diversity that different models produce, with specific emphasis and different aspects of the features they want to model (Zhou et al., 2020a). In the last few decades, many approaches have been developed to forecast streamflow, including physically based and data-driven models (Tikhamarine et al., 2020; Zuo et al., 2020). Although physically based models can help understand underlying physical processes, they usually require a large amount of input information, such as meteorological and geographic data as well as soil and land use characteristics (Guo et al., 2018, 2020a). Different from physically based models, data-driven models based on statistical modeling have attracted significant interest due to their simplicity and satisfactory forecast results with low information requirements (Al-Sudani et al., 2019; Mehdizadeh et al., 2019; Osman et al., 2020). Artificial intelligence (AI)-based approaches, i.e., machine learning (ML) methods, belong to the latter group. Widely used ML approaches include artificial neural networks (ANNs) and least-squares support vector machines (LSSVMs) (Ghumman et al., 2018; Kisi et al., 2019; Meng et al., 2019; Adnan et al., 2020; Ali and Shahbaz, 2020). Such models have been proven to be efficient tools to model qualitative and quantitative hydrological variables and deal with nonlinear features in streamflow. In recent years, the booming development of deep learning technology has brought many new approaches, such as recurrent neural networks (RNNs) (Elman, 1990), one of the most popular neural networks in the deep learning field. RNNs can preserve and remember the short-past and long-past information and thus are preferred for a complex and highly nonlinear timing problem. Long short-term memory (LSTM) (Hochreiter and Schmidhuber, 1997) and gated recurrent unit (GRU) models (Cho et al., 2014) are two different versions of RNNs. LSTM and GRU networks have been successfully applied in many fields (Greff et al., 2017; Zhang et al., 2018; Jung et al., 2020; Shahid et al., 2020; Ayzel and Heistermann, 2021), and they are demonstrated to generate comparable performances, But GRU has a more straightforward structure and a higher operation speed than LSTM. Recently, many applications that assessed them together are also found in the hydrological field (Gao et al., 2020; Muhammad et al., 2019).

While a considerable research effort has been made to evaluate and improve the quality of streamflow forecasts (Gibbs et al., 2018; Nanda et al., 2019; Sharma et al., 2019; Van Osnabrugge et al., 2019; Feng et al., 2020; Pechlivanidis et al., 2020), how forecasts impact decision-making in the real-time reservoir operations has also gradually gained researchers' attention (Goddard et al., 2010; Shamir, 2017; Anghileri et al., 2019; Alexander et al., 2020; Hadi et al., 2020), e.g., do high-quality forecasts mean improved decision? Traditionally, a skillful forecast is vital for the reliability of the forecasts and is essential to promote the use of forecasts in real-world applications by decision makers. In fact, forecast value is expected to increase with forecast quality, but it may also vary based on other factors such as reservoir capacity and operating objectives (Anghileri et al., 2016). Some studies have even disproved the intuitive assumption that higher forecast performance always leads to better operation decisions, for example, in agricultural water management (Chiew et al., 2003) and water resources allocation (Turner et al., 2017). Therefore, when forecasts are used to support reservoir operation, it should be assessed in which conditions they can help make better decisions. Moreover, forecast uncertainty and error generally grow with the increase of the forecast horizon (Maurer and Lettenmaier, 2004; Denaro et al., 2017; Zhao et al., 2019). A decision maker may doubt whether longer forecast lead times provide sufficient information for a decision purpose or not. There is often a mismatch between the information needed for reservoir operations and the skillful lead time of the reservoir inflow forecast (Anghileri et al., 2016). It is crucial to demonstrate the applicability and effectiveness of the forecast horizon in a forecast-based reservoir operation system (Xu et al., 2014). Overall, there is a continuous need for in-depth study to conduct posterior evaluations of forecasts with different forecast lead times and obtain the efficient forecast horizon for water allocation.

A decision maker must allocate limited water to different water use sectors considering the conflicting objectives (e.g., benefits and costs) and multiple uncertainties (e.g., forecast uncertainty) in a forecast-based reservoir operation system. Multi-objective programming (MOP) is a valuable tool for helping decision makers facilitate decision-making with multiple conflicting objectives (Fang et al., 2018b; Guo et al., 2020c), which can offer feasible methods for generating compromise decision alternatives. Some MOP approaches have been widely developed to tackle the uncertainty associated with the decision-making processes, such as multi-objective fuzzy programming (Zimmermann, 1978; Pishvaee and Razmi, 2012; Ren et al., 2017) and multi-objective stochastic programming (Xu et al., 2014, 2020; Zhang et al., 2020). These approaches generally convert the multi-objective functions into a single-objective deterministic problem through a fuzzy programming method or a constraint operator. They can effectively deal with the uncertainties between objectives and/or constraints by integrating the decision makers aspiration levels. However, they may encounter difficulties due to the need for predetermined individual preferences or reasonable bounds for all objectives. In comparison, multi-objective robust decision-making (MORDM) is an effective way to handle such difficulties (Kasprzyk et al., 2013; Zeff et al., 2014; Yan et al., 2017; Hadjimichael et al., 2020). It can generate many alternative solutions (Pareto solutions) that do not require assumptions about decision makers' preferences and enhance the robustness of the optimization process. Besides, MORDM, by parameterizing the decision space, can avoid the curse of dimensionality in some MOP approaches, simplify computational complexity, and reduce the running time (Giuliani et al., 2016; Salazar et al., 2017).

Figure 1 Framework of the AI-based management methodology.

In summary, there are still several challenges in forecast-informed reservoir optimization. To address these challenges, the specific research questions of this study are as follows:

1. Can GRU achieve the same accuracy in the streamflow forecast compared to LSTM with fewer parameters and more straightforward structures?

2. In which conditions can an improvement in forecast skill be translated into improved reservoir operation optimization?

3. How do such short-term inflow forecasts with different forecast horizons be used to optimize the multi-reservoir system to impact operation results?

To answer the questions mentioned above, we build an AI-based management framework, which integrates multi-step streamflow forecasts and multi-reservoir operation optimization. We strive to (1) simulate inflow using LSTM, GRU, and LSSVM and verify their effectiveness on short-term deterministic streamflow forecasts; (2) generate stochastic inflow scenarios using BMA for refining uncertainty characterization; (3) develop the parameterized MORDM framework for a multi-reservoir operation system and inform decision-making by assessing the value, that is, the operation benefit gain or the induced cost of forecasts with a particular lead time. As a case study, including one recipient reservoir storing water from the continental diversion project and 24 supply reservoirs storing water from local rainfall, 25 reservoirs supplying water for four water plants in the Zhoushan Islands, China, are chosen to assess the performance of the AI-based forecast and the forecast-informed operation.

2 Methodology

The experimental approach followed in the study is shown in Fig. 1 and described in the following sections.

2.1 Machine learning (ML) methods

This section gives a brief introduction to long short-term memory (LSTM), gated recurrent unit (GRU), and least-squares support vector machine (LSSVM) models. In this study, the mapping function between the forecasted streamflow Q_t and hydrological variables x_t can be represented by f(⋅). In LSTM and GRU, $Q_t=f(x_t,h_{t-\mathrm{1}})\left(h_{t-\mathrm{1}}\right)$ denotes the last hidden cell state and the initial state of h_t is h₀=0, while in GWO-LSSVM, Q_t=f(x_t).

Figure 2 Structure of an (a) LSTM and (b) GRU cell.

2.1.1 Long short-term memory (LSTM)

The LSTM network is one of the recurrent neural networks (RNNs) developed by Hochreiter and Schmidhuber (1997), and the basic structure of an LSTM cell is illustrated in Fig. 2a. It is an improved RNN aiming to solve problems such as gradients in long-term memory and backpropagation. The LSTM cell has three gates maintaining and adjusting its cell state and hidden state, including the forget gate, input gate, and output gate. The forget gate determines what information should be thrown away from the cell state. The input gate decides which information is used to update the cell state. The output gate controls which information stored in the current cell state flows into the new hidden state. In Fig. 2a, the state (c_t) and the hidden state (h_t) of the LSTM cell are updated as follows (Hochreiter and Schmidhuber, 1997):

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathrm{forget}\mathrm{gate},g_t=\mathit{\sigma }\left(W_f{x}_t+U_f{h}_{t-\mathrm{1}}+b_f\right);\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{input}\mathrm{gate},i_t=\mathit{\sigma }\left(W_i{x}_t+U_i{h}_{t-\mathrm{1}}+b_i\right);\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{potential}\mathrm{cell}\mathrm{state},{\stackrel{\mathrm{̃}}{c}}_t=\tanh \left(W_{\stackrel{\mathrm{̃}}{c}}{x}_t+U_{\stackrel{\mathrm{̃}}{c}}{h}_{t-\mathrm{1}}+b_{\stackrel{\mathrm{̃}}{c}}\right);\text{(4)}& {\displaystyle }& {\displaystyle }\mathrm{cell}\mathrm{state},c_t=f_t\odot c_{t-\mathrm{1}}+i_t\odot {\stackrel{\mathrm{̃}}{c}}_t;\text{(5)}& {\displaystyle }& {\displaystyle }\mathrm{output}\mathrm{state},o_t=\mathit{\sigma }\left(W_o{x}_t+U_o{h}_{t-\mathrm{1}}+b_o\right);\text{(6)}& {\displaystyle }& {\displaystyle }\mathrm{hidden}\mathrm{state},h_t=o_t\tanh \odot c_t;\end{array}$$

where g_t, c_t, o_t, and ${\stackrel{\mathrm{̃}}{c}}_t$ represent the forget gate, input gate, output gate, and potential cell state, respectively. ⊙ denotes the element-wise multiplication of vectors, and tanh(⋅) is the hyperbolic tangent. x_t represents the current input vector, h_t−1 denotes the last hidden cell state, and the initial state of h_t is h₀=0. σ(⋅) represents the logistic sigmoid function. [W_f, W_i, W_o, $W_{\stackrel{\mathrm{̃}}{c}}$], [U_f, U_i, U_o, $U_{\stackrel{\mathrm{̃}}{c}}$], and [b_f, b_i, b_o, $b_{\stackrel{\mathrm{̃}}{c}}$] represent the input weight matrix, recurrent weight matrix, and bias vectors for the forget, input–output, and potential cell gates, respectively.

2.1.2 Gated recurrent unit (GRU)

GRU networks were proposed as a modification of LSTM networks with a more straightforward structure (Cho et al., 2014). The specific structure of the GRU cell is shown in Fig. 2b. Compared with LSTM, GRU only has two control gates, including a reset gate and an update gate. The update gate is applied to control how much information of the previous step is brought into the current step, while the reset gate is used to control the degree of ignoring the information of the previous state. In this way, GRU is superior to LSTM in terms of computer modeling time and parameter updates. In Fig. 2b, the state (c_t) and the hidden state (h_t) of the GRU cell are updated as follows (Cho et al., 2014):

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }\mathrm{reset}\mathrm{state},r_t=\mathit{\sigma }\left(W_r{x}_t+U_r{h}_{t-\mathrm{1}}+b_r\right);\text{(8)}& {\displaystyle }& {\displaystyle }\mathrm{update}\mathrm{state},z_t=\mathit{\sigma }\left(W_z{x}_t+U_z{h}_{t-\mathrm{1}}+b_z\right);\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{potential}\mathrm{cell}\mathrm{state},{\stackrel{\mathrm{̃}}{c}}_t=\tanh \left(W_{\stackrel{\mathrm{̃}}{c}}{x}_t+U_{\stackrel{\mathrm{̃}}{c}}\left(r_t\odot h_{t-\mathrm{1}}\right)+b_{\stackrel{\mathrm{̃}}{c}}\right);\text{(10)}& {\displaystyle }& {\displaystyle }\mathrm{cell}\mathrm{state},c_t=\left(\mathrm{1}-z_t\right)\odot c_{t-\mathrm{1}}+z_t\odot {\stackrel{\mathrm{̃}}{c}}_t;\text{(11)}& {\displaystyle }& {\displaystyle }\mathrm{hidden}\mathrm{state},h_t=c_t;\end{array}$$

where r_t, z_t, and ${\stackrel{\mathrm{̃}}{c}}_t$ represent the reset, update, and potential cell state, respectively. ⊙ denotes the element-wise multiplication of vectors, and tanh(⋅) is the hyperbolic tangent. x_t represents the input vectors, h_t−1 denotes the last hidden cell state, and the initial state of h_t is h₀=0. σ(⋅) represents the logistic sigmoid function. [W_r, W_z, $W_{\stackrel{\mathrm{̃}}{c}}$], [U_r, U_z, $U_{\stackrel{\mathrm{̃}}{c}}$], and [b_r, b_z, $b_{\stackrel{\mathrm{̃}}{c}}$] represent the input weight matrix, recurrent weight matrix, and bias vectors for the reset, update, and potential cell gates, respectively.

2.1.3 Least-squares support vector machine with grey wolf optimizer (GWO-LSSVM)

LSSVM is a modified version of SVM, proposed by Suykens and Vandewalle (1999) to reduce the computational time of SVM. SVM uses the quadratic program to formulate the training process of the modeling procedure, while LSSVM aims to employ the least-squares loss functions. The LSSVM nonlinear function is expressed as (Suykens et al., 2002)

$$\begin{array}{}\text{(12)}& f\left(x\right)=w^T\mathit{\phi }\left(x\right)+b,\end{array}$$

where φ(⋅) is the mapping function that maps the input x into a d-dimensional feature vector, w is a weight vector, and b represents bias. In LSSVM, a minimum objective function is proposed to estimate ω and b (Suykens et al., 2002).

$$\begin{array}{}\text{(13)}& \min J(w,e)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^{T}w+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\gamma }\sum _{i=\mathrm{1}}^N{e}_i^{\mathrm{2}},\end{array}$$

that has the following constraints (Suykens et al., 2002):

$$\begin{array}{}\text{(14)}& y_i=w^T\mathit{\phi }\left(x_i\right)+b+e_i,\end{array}$$

where e is the error variable, and γ is the regulative constant. The objective function can be obtained to solve the optimization problems in Eq. (15) by introducing the Lagrange multipliers α and transferring the constraint problem into an unconstrained one (Suykens et al., 2002):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}L(w,b,e,\mathit{\alpha })& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{w}^{T}w+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathit{\gamma }\sum _{i=\mathrm{1}}^N{e}_i^{\mathrm{2}}\\ & -\sum _{i=\mathrm{1}}^N{\mathit{\alpha }}_i\left(w^T\mathit{\phi }\left(x_i\right)+b+e_i-y_i\right).\end{array}\end{array}$$

By finding the partial derivative of Eq. (16) with respect to w, b, α_i, and e_i, the following equation can be derived:

$$\begin{array}{}\text{(16)}& y=\sum _{i=\mathrm{1}}^N{\mathit{\alpha }}_i\left(\mathit{\phi }(x{)}^T\mathit{\phi }\left(x_i\right)\right)+b=\sum _{i=\mathrm{1}}^N{\mathit{\alpha }}_{i}K\left(x,x_i\right)+b,\end{array}$$

where K(x,x_i) is the kernel function. Many kernel functions such as linear, polynomial, radial basis, and sigmoidal have been proposed for LSSVM (Bemani et al., 2020). We adopt the most widely used kernel function, radial basis function (RBF), in this study. The RBF is expressed as

$$\begin{array}{}\text{(17)}& K\left(x,x_i\right)=\exp \left(-{∥x-x_i∥}^{\mathrm{2}}/\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\right),\end{array}$$

where σ² is the kernel parameter. In this study, the parameters γ and σ were optimized using the grey wolf optimizer (GWO). Please see more details on GWO in Guo et al. (2020d).

2.2 Bayesian model averaging (BMA)

Generally, it is difficult to determine which model is the best one, leading to model uncertainty. BMA is proposed to solve the uncertainty of the models through averaged estimations from individual models (Liu and Merwade, 2019; Samadi et al., 2020). The weight for each model is based on the simulated decision probability density function, i.e., the posterior probability of the model. Suppose Q is the unknown quantity we want to predict, given a subset of model forecasts $f=\mathit{\{}{f}_{\mathrm{1}},f_{\mathrm{2}}$, …, f_K} (k=1, 2, …, K, where K is the number of individual model) and the observed data D, the posterior distribution of Q can be calculated as (Hoeting et al., 1999)

$$\begin{array}{}\text{(18)}& \begin{array}{rl}p\left(Q\right|D)& =\sum _{i=\mathrm{1}}^{K}p\left(f_k|D\right)\cdot p_k\left(Q|f_k,D\right)\\ & =\sum _{i=\mathrm{1}}^K{w}_k\cdot p_k\left(Q|f_k,D\right),\end{array}\end{array}$$

where $p_k\left(Q\right|f_k,D)$ is the posterior distribution of Q given the model forecast f_k and the observed data D, and p(f_k|D) is the posterior probability. In this case, posterior probabilities are the weighting factor for each model, and $\sum _{k=\mathrm{1}}^K{w}_k=\mathrm{1}$. The posterior mean (E) and variance (V) of Q are as follows (Hoeting et al., 1999):

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }E\left[Q\right|D]=\sum _{k=\mathrm{1}}^K{w}_k\cdot E\left[p_k\left(Q|f_k,D\right)\right]=\sum _{k=\mathrm{1}}^K{w}_k{f}_k,\text{(20)}& {\displaystyle }& {\displaystyle }V\left[Q\right|D]=\sum _{k=\mathrm{1}}^K{w}_k\cdot {\left[f_k-\sum _{k=\mathrm{1}}^K{w}_k{f}_k\right]}^{\mathrm{2}}+\sum _{k=\mathrm{1}}^K{w}_k{\mathit{\sigma }}_k^{\mathrm{2}},\end{array}$$

where w_k and ${\mathit{\sigma }}_k^{\mathrm{2}}$ are the posterior mean (weight) and variance of the kth forecast model. In this study, a log-likelihood function is maximized to estimate the parameters (weight w_k and variance ${\mathit{\sigma }}_k^{\mathrm{2}}$) as shown in Eq. (21).

$$\begin{array}{}\text{(21)}& l\left(\mathit{\theta }\right)=\log \left(\sum _{k=\mathrm{1}}^K\left(w_k\cdot g\left(Q|f_k^t,{\mathit{\sigma }}_k^{\mathrm{2}}\right)\right)\right),\end{array}$$

where θ is the vector of parameters {w_k, ${\mathit{\sigma }}_k^{\mathrm{2}}$, k=1, 2, …, K}. $g\left(Q\right|f_k^t$, ${\mathit{\sigma }}_k^{\mathrm{2}})$ is the Gaussian distribution function, where w_k is the weight and ${\mathit{\sigma }}_k^{\mathrm{2}}$ is the variance.

The expectation–maximization (EM) algorithm (Lee et al., 2020) is used to find out the maximum likelihood with a termination criterion (early stopping or a maximal iteration). As the EM proceeds, the parameters of weight w_k and variance ${\mathit{\sigma }}_k^{\mathrm{2}}$ are updated as follows.

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{w}_k^{\left(\mathrm{Iter}\right)}={\displaystyle \frac{\mathrm{1}}{NT}}\left(\sum _{t=\mathrm{1}}^{NT}{z}_k^{t^{\left(\mathrm{Iter}\right)}}\right),\text{(23)}& {\displaystyle }& {\displaystyle }{\mathit{\sigma }}_k^{{\mathrm{2}}^{\left(\mathrm{Iter}\right)}}={\displaystyle \frac{\sum _{t=\mathrm{1}}^{NT}{z}_k^{t^{\mathrm{Iter}}}\cdot {\left(Y^t-f_k^t\right)}^{\mathrm{2}}}{\sum _{t=\mathrm{1}}^{NT}{z}_k^{t^{\mathrm{Iter}}}}},\text{(24)}& {\displaystyle }& {\displaystyle }{z}_k^{t^{\left(\mathrm{Iter}\right)}}={\displaystyle \frac{g\left(Q|f_k^t,{\mathit{\sigma }}_k^{{\mathrm{2}}^{(\mathrm{Iter}-\mathrm{1})}}\right)}{\sum _{k=\mathrm{1}}^{K}g\left(Q|f_k^t,{\mathit{\sigma }}_k^{{\mathrm{2}}^{(\mathrm{Iter}-\mathrm{1})}}\right)}},\text{(25)}& {\displaystyle }& {\displaystyle }l(\mathit{\theta }{)}^{\left(\mathrm{Iter}\right)}=\sum _{t=\mathrm{1}}^{NT}\log \left(\sum _{k=\mathrm{1}}^K\left(w_k^{\left(\mathrm{Iter}\right)}\cdot g\left(Q|f_k^t,{\mathit{\sigma }}_k^{{\mathrm{2}}^{\left(\mathrm{Iter}\right)}}\right)\right)\right),\end{array}$$

where “Iter” is the number of iterations. NT is the length of calibration periods. Y^t and $f_k^t$ are the observed and forecast streamflow at the tth time step, respectively (m³ s⁻¹), and $z_k^{t^{\left(\mathrm{Iter}\right)}}$ is the latent variable for the kth model at the tth time step in the Iter iteration. Then we use the Monte Carlo simulation method to generate BMA ensemble forecasts. Assume M is the number of Monte Carlo simulations, and the procedure is described below (Zhou et al., 2020a).

• a.

  Set the initial cumulative weight $w_{\mathrm{0}}^{*}=\mathrm{0}$, and calculate the cumulative weight $w_k^{*}=w_{k-\mathrm{1}}^{*}+w_k$ for k=1, 2, …, K. Create a random variable u between 0 and 1. If $w_{k-\mathrm{1}}^{*}\le u\le w_k^{*}$, the kth forecast model would be used as the target forecast.

• b.

  Generate a realization of the forecasts Q_t using the Gaussian distribution function $g\left(Q_t\right|f_k^t$, ${\mathit{\sigma }}_k^{\mathrm{2}})$. In such a way, there are a set of alternative forecasts to be chosen from as the final forecast.

• c.

  Repeat steps (a) and (b) M times, and obtain a set of streamflow series $Q_{T\times M}=\left[\begin{array}{cccc}{Q}_{\mathrm{1},\mathrm{1}}& Q_{\mathrm{2},\mathrm{1}}& \mathrm{\dots }& Q_{T,\mathrm{1}}\\ Q_{\mathrm{1},\mathrm{2}}& Q_{\mathrm{2},\mathrm{2}}& \mathrm{\dots }& Q_{T,\mathrm{2}}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ Q_{\mathrm{1},M}& Q_{\mathrm{1},M}& \mathrm{\dots }& Q_{T,M}\end{array}\right]$. Furthermore, 90 % confidence intervals between the 5 % and 95 % quantities are employed to represent the uncertainty of BMA ensemble forecasts.

2.3 Forecast performance measures

Three performance indicators are applied to assess the deterministic forecast performance of the three data-process models. They are the Nash–Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970), the root mean square error (RMSE) (Karunanithi et al., 1994), and the mean absolute error (MAE) (Legates and McCabe, 1999). They are expressed as below:

$$\begin{array}{}\text{(26)}& {\displaystyle }& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{\sum _{t=\mathrm{1}}^T{\left(Q_{m,t}-Q_{o,t}\right)}^{\mathrm{2}}}{\sum _{t=\mathrm{1}}^T{\left(Q_{o,t}-{\stackrel{\mathrm{‾}}{Q}}_o\right)}^{\mathrm{2}}}},\text{(27)}& {\displaystyle }& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T{\left(Q_{m,t}-Q_{o,t}\right)}^{\mathrm{2}}},\text{(28)}& {\displaystyle }& {\displaystyle }\mathrm{MAE}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T\left|Q_{m,t}-Q_{o,t}\right|,\end{array}$$

where T is the number of samples, Q_m,t is the forecasted reservoir inflow (m³ s⁻¹), Q_o,t is the observed inflow (m³ s⁻¹), and ${\stackrel{\mathrm{‾}}{Q}}_o$ is the average of the observed inflow (m³ s⁻¹). The NSE can be used to evaluate the stability of the forecasted value. In contrast, the RMSE and MAE are used to characterize the overall forecast accuracy. The NSE value is (−∞, 1], while MAE and RMSE values are (0, +∞). Generally, models with larger values of NSE or smaller values of RMSE and MAE provide better forecasting accuracy.

In addition, two performance indicators are used to evaluate the performance of ensemble forecast models, i.e., the containing ratio (CR), and average deviation amplitude (D), which were adopted for assessing the goodness of the prediction bounds (Xiong et al., 2009).

$$\begin{array}{}\text{(29)}& {\displaystyle }\begin{array}{rl}& \mathrm{CR}={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T{N}_t\times \mathrm{100}\mathit{\%}\\ & N_t=\left\{\begin{array}{ll}\mathrm{1}& \mathrm{if}{\widehat{Q}}_{l,t}\le Q_{o,t}\le {\widehat{Q}}_{u,t}\\ \mathrm{0}& \mathrm{else}\end{array}\right.,\end{array}\text{(30)}& {\displaystyle }D={\displaystyle \frac{\mathrm{1}}{T}}\sum _{t=\mathrm{1}}^T\left|{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\widehat{Q}}_{l,t}+{\widehat{Q}}_{u,t}\right)-Q_{o,t}\right|,\end{array}$$

where ${\widehat{Q}}_{l,t}$ and ${\widehat{Q}}_{u,t}$ represent the lower and upper prediction bounds of streamflow (m³ s⁻¹), respectively. Clearly, models with higher CR values but lower D values would produce better performance.

2.4 Parameterized multi-objective robust decision-making (MORDM)

This study proposes a parameterized multi-objective robust decision-making approach to design operating policies for the multi-objective reservoir operations by combining direct policy search (DPS) and multi-objective robust decision-making (MORDM). In the parameterized MORDM, instead of using the volumes of water to be allocated as the decision variables, we prescribe decisions approximated as nonlinear functions conditioned on system state variables (e.g., forebay water level observed or predicted inflows and precipitation) (Giuliani et al., 2016; Quinn et al., 2017a, b; Salazar et al., 2017). The nonlinear functions can be realized by the DPS approach. DPS is based on the parameterization of the operating policy p_θ and the exploration of the parameter space Θ to find a parameterized policy that optimizes the expected function, i.e.,

$$\begin{array}{}\text{(31)}& p_{\mathit{\theta }}^{*}={\mathrm{argmin}}_{p_{\mathit{\theta }}}{\left(J_{\mathrm{1}},J_{\mathrm{2}},\mathrm{\dots },J_M\right)}_{p_{\mathit{\theta }}}s\cdot t\cdot \mathit{\theta }\in \mathrm{\Theta },\end{array}$$

where J₁, J₂, …, J_M are the objective functions, and M is the number of objectives. $p_{\mathit{\theta }}^{*}$ is the corresponding optimal policy with parameters θ^*. Different DPS approaches have been proposed, where two nonlinear approximating networks, namely artificial neural networks (ANNs) and radial basis functions (RBFs), have become widely adopted as universal approximators in many applications (Deisenroth et al., 2013). In particular, we parameterize the operating policy as RBFs because they have been demonstrated to be effective in solving multi-objective water resources management problems (Giuliani et al., 2014, 2015), and the kth decision variable in the vector u_t (with k=1, …, K) is defined as

$$\begin{array}{}\text{(32)}& u_t^k=\sum _{i=\mathrm{1}}^N{\mathit{\omega }}_{i,k}{\mathit{\phi }}_{i,k}\left({\mathrm{\Gamma }}_t\right),\end{array}$$

where N is the number of RBFs φ(⋅), Γ_t is the policy input vectors at the tth time step including exogenous information (e.g., forebay water level observed or predicted inflows and precipitation), and ω_i,k is the weight of the ith RBF, $\sum _{i=\mathrm{1}}^N{\mathit{\omega }}_{i,k}=\mathrm{1}{\mathit{\omega }}_{i,k}>\mathrm{0}$. The single RBF is defined as follows:

$$\begin{array}{}\text{(33)}& {\mathit{\phi }}_{i,k}\left({\mathrm{\Gamma }}_t\right)=\exp \left[-\sum _{j=\mathrm{1}}^M{\displaystyle \frac{{\left[{\left({\mathrm{\Gamma }}_t\right)}_j-c_{j,i}\right]}^{\mathrm{2}}}{b_{j,i}^{\mathrm{2}}}}\right],\end{array}$$

where (Γ_t)_j is the jth policy input at the tth time step, and M denotes the number of policy input vectors Γ_t, j=1, 2, …, M. c_i and b_i are the M-dimensional center and radius vectors of the ith RBF, respectively. The centers of the RBF must lie within the bounded input space (Yang et al., 2017). The parameter vector θ is defined as $\mathit{\theta }=[c_{i,j,k},b_{i,j,k},{\mathit{\omega }}_{i,j,k}]$, where the number of θ is $n_{\mathit{\theta }}=N\times K\times (\mathrm{2}\times M+\mathrm{1})$. In general, when DPS problems involve multiple objectives, they can be coupled with truly multi-objective optimization methods, such as multi-objective evolutionary algorithms (MOEAs), which allow for an approximation of the Pareto front in a single run of the algorithm.

Figure 3 Schematization of the parameterized MORDM methods.

In our study, the parameterized MORDM approach will be coupled with a rolling horizon scheme over a 1-year period to solve the multi-objective reservoir operation problem. Given the lead time of 7 d (the forecast horizon is equal to the operation horizon) as an example, it is operated following two steps: the optimization model is first operated daily over a 7 d horizon using the parameterized MORDM; after implementing current water allocation decisions, the status, inflow, and other information of reservoirs update as time evolves, and then the remainder is subsequently operated. The two steps are repeated until the process (1-year period) is completed. In each operating horizon, the main steps of the parameterized MORDM are described below and presented in Fig. 3.

1. Problems are formulated, including the performance measures and constraints.

2. Generate alternative parameterized policies subject to all the constraints, and the objectives are evaluated over stochastic inflows with the following procedures (Giuliani et al., 2016):

   • a.

     The operating policies are parameterized using RBFs.

   • b.

     Run a system simulation from t=1, 2, … 7 d upon each individual parameterized policy p_θ for each inflow series and obtain the system trajectories.

   • c.

     Compute the parameterized policies performance in terms of the operating objectives as a function of system trajectories.

3. Recompute the parameterized policy performance with robust criteria, for instance, the principle of insufficient reason, minimax, and minimax regret (Guo et al., 2020b). Among them, the principle of insufficient reason transforming the problem under uncertainty into a decision-making problem under risk has been used in water resources problems (Giuliani and Castelletti, 2016). The principle of insufficient reason suggests that in the absence of knowledge on the probabilities associated with the different states, the decision could be taken by assigning equal probability to all the states (i.e., $P_j=\mathrm{1}/n$). The robust parameterized policy performance can be expressed as

   $$\begin{array}{}\text{(34)}& \min \left({\displaystyle \frac{\mathrm{1}}{n}}\sum _{j=\mathrm{1}}^n\mathrm{Obj}\left(p_{\mathit{\theta }},s_j\right)\right),\end{array}$$

   where Obj(p_θ,s_j) is the performance function using parameterized policy p_θ upon jth streamflow series, s_j denotes the scenario of the jth streamflow series, and n is the number of stochastic streamflow series.

4. Optimize the parameterized policies using multi-objective evolutionary algorithms (MOEAs) based on the robust performance objectives. Repeat steps (2)–(4) until the times of population iteration are reached, and export the optimal Pareto solutions. In this study, the optimization is solved by applying NSGA-II to search the space of decision variables and identify the trajectories.

It should be noted that the parameterized MORDM in this study aims to solve optimization problems under uncertainty, and thereby one streamflow series needs to be repeated multiple times.

Figure 4 Location of the Zhoushan Islands.

3 Case study

3.1 Study area and data

The Zhoushan Islands are located in the northeast of Zhejiang Province, China, with a total area of 22 000 km² and 1390 islands (Fig. 4). The climate is governed by monsoon-influenced subtropical marine weather systems, and the annual mean temperature and precipitation are 17 ^∘C and 1300 mm, respectively. There are no large rivers in the islands, and the insufficient freshwater resources severely limit the development of industry and population in the Zhoushan Islands. Recently, a continental diversion project transferring water from the city of Ningbo to Zhoushan Islands has been treated as an effective solution to overcome the water scarcity problem partially. The transferred water is stored in Huangjinwan Reservoir and then operated together with the limited freshwater resources in the remaining 24 reservoirs to supply water to four water plants, i.e., Daobei, Hongqiao, Lincheng, and Pingyangpu. Data for this study include historical inflow and state of reservoirs, water demand of water plants, and climate forcing data over 2002–2008. The climate data, including daily precipitation and evaporation, are observed at one meteorological station and three rainfall stations. The characteristics of the reservoirs are listed in Table 1.

3.2 Problem formulation

Figure 5 shows the simplified schematic diagram of the water supply system in Zhoushan Islands, including reservoirs, pumping stations, pipelines, and water plants. The pipeline arrow indicates the direction of the water flow. It covers the processes associated with water abstraction from resources, water distribution through the network involving the use of pumping stations and pipelines, and main activities relevant to water flow. In this study, water resources include local surface water and imported water. The surface water is the water stored in local reservoirs (a number of 24 reservoirs), while the imported water is the water transferred from the city of Ningbo (stored in Huangjinwang Reservoir). The imported water is transferred from the city of Ningbo to Zhoushan Islands through Lixidu and Lanshan pumping stations. End users within the water supply system are generally divided into the household, industry, agriculture, and environmental use. This study mainly considers household and industry use, which water plants can supply. The agriculture and environmental use are satisfied through operating the reservoir storage above a specific value. That is to say, the main goal of the water allocation plan is to ensure sufficient water flows into the four plants in Zhoushan Islands. They are Daobei, Lincheng, Hongqiao, and Pingyangpu plants, respectively. Releases from the reservoirs (Huangjinwan Reservoir and the remaining 24 local reservoirs) must meet the requirements of water plants. As observed in Fig. 5, the reservoirs supplying plants can be divided into two categories. Some reservoirs can directly release water into the plants or reservoirs, including Longtan, Ludong, Shatianao, Nanao, Chenao, Cengang, Tuanjie, and Changchunling reservoirs. In contrast, the other reservoirs can only release water into the plants or reservoirs using pumping stations. In such a way, the pumping flow can be obtained by summing reservoir releases through the corresponding pumping station, using the following equation.

$$\begin{array}{}\text{(35)}& Q_{t,j}^p=\sum _{n=\mathrm{1}}^{N_{\mathrm{1}}}{Q}_{t,n}^r,\end{array}$$

where $Q_{t,j}^p$ denotes the jth pumping flow at the tth time step (m³ s⁻¹), $Q_{t,n}^r$ denotes the release of the nth reservoir at the tth time step (m³ s⁻¹), and N₁ is the number of reservoirs pumped by the jth pumping station.

Figure 5 Schematic diagram of the Zhoushan Islands.

It can be noted in Fig. 5 that there are no specific hydraulic connections between most of the reservoirs, while Chahe, Hongwei, Chengbei, and Xiamen reservoirs can release water into Hongqiao Reservoir (the largest reservoir in Zhoushan Islands). With a water plant as a center, the whole islands are divided into four districts, i.e., Daobei, Lincheng, Hongqiao, and Dongbu. The dashed line represents the district boundary. Each district includes a water plant, several pumping stations, and reservoirs to supply water for the water plant. The hydraulic connection between such a water plant and corresponding pumping stations and reservoirs can be expressed as

$$\begin{array}{}\text{(36)}& W_t^s=\sum _{j=\mathrm{1}}^J{Q}_{t,j}^p\mathrm{\Delta }t+\sum _{n=\mathrm{1}}^{N_{\mathrm{2}}}{Q}_{t,n}^r\mathrm{\Delta }t,\end{array}$$

where $W_t^s$ is the amount of water supply for a water plant at the tth time step (m³), J is the number of pumping stations flowing into the water plant, and N₂ is the number of reservoirs directly releasing into the water plant.

Table 1 Reservoir characteristics in the Zhoushan Islands.

In Fig. 5, every two system elements are connected by the pipelines, e.g., reservoir and reservoir, reservoir and pumping station, and pumping station and water plant. In some cases, more than one reservoir or pumping station share one pipeline, leading to competition on channel flow. However, the multi-objective optimization problem is operated on a daily time step in our study, and we assume reservoir releases or pumping station flows to the water plant without considering the channel flow limitation, thereby regardless of the specific hydrologic connections between channels or pipelines.

Three objectives are identified to evaluate the performance of the strategies. The conflicting objectives are to minimize the water deficiency ratio of the Daobei Plant, minimize the water deficiency ratio of the remaining three plants (Hongqiao, Lincheng, and Pingyangpu), and maximize the net benefits. The three plants can feed each other and thus are considered together in our study. A decision maker would consider a different suite of costs depending on whether an existing system is being managed or a completely new system is being designed. As water supply occurs in an existing system, costs considered in this study are the operating costs. These objective functions are given as follows:

$$\begin{array}{}\text{(37)}& {\displaystyle }\begin{array}{rl}\mathrm{Min}{\mathrm{Obj}}_{\mathrm{1}}\left(x\right)=& \left(\sum _{t=\mathrm{1}}^T{W}_t^{n,\mathrm{db}}-\sum _{t=\mathrm{1}}^T{W}_t^{s,\mathrm{db}}\right)/\\ & \sum _{t=\mathrm{1}}^T{W}_t^{n,\mathrm{db}}\times \mathrm{100}\mathit{\%},\end{array}\text{(38)}& {\displaystyle }\begin{array}{rl}\mathrm{Min}{\mathrm{Obj}}_{\mathrm{2}}\left(x\right)=& \left(\sum _{k=\mathrm{1}}^{\mathrm{3}}\sum _{t=\mathrm{1}}^T{W}_{t,k}^{n,\mathrm{th}}-\sum _{k=\mathrm{1}}^{\mathrm{3}}\sum _{t=\mathrm{1}}^T{W}_{t,k}^{s,\mathrm{th}}\right)/\\ & \sum _{i=\mathrm{1}}^{\mathrm{3}}\sum _{t=\mathrm{1}}^T{W}_{t,k}^{n,\mathrm{th}}\times \mathrm{100}\mathit{\%},\end{array}\text{(39)}& {\displaystyle }\mathrm{Min}{\mathrm{Obj}}_{\mathrm{3}}\left(x\right)=\left(M_{\mathrm{c}}^{\mathrm{lw}}+M_{\mathrm{c}}^{\mathrm{iw}}\right)-M_{\mathrm{r}},\end{array}$$

where Obj₁ and Obj₂ are the water deficiency ratio of Daobei Plant and the sum of the remaining three plants, respectively (%); Obj₃ denotes the net operating costs (RMB); $W_t^{s,\mathrm{db}}$ and $W_t^{n,\mathrm{db}}$ are the amount of water supply and demand for Daobei Plant at the tth time step, respectively (m³); $W_{t,k}^{n,\mathrm{th}}$ and $W_{t,k}^{n,\mathrm{th}}$ are the amount of water supply and demand for the kth plant (one of the remaining three plants) at the tth time step, respectively (m³); $M_{\mathrm{c}}^{\mathrm{lw}}$ and $M_{\mathrm{c}}^{\mathrm{iw}}$ are the operating costs for water supply using local reservoir water and imported water, respectively (RMB); and M_r is the revenue (RMB). The costs and revenue can be obtained according to the following:

1. operating costs for water supply using local reservoir water ($M_{\mathrm{c}}^{\mathrm{lw}}$, RMB),

   $$\begin{array}{}\text{(40)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c}}^{\mathrm{lw}}=M_{\mathrm{c},\mathrm{1}}^{\mathrm{lw}}+M_{\mathrm{c},\mathrm{2}}^{\mathrm{lw}}+M_{\mathrm{c},\mathrm{3}}^{\mathrm{lw}}\text{(41)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{1}}^{\mathrm{lw}}=c_{\mathrm{1}}^{\mathrm{lw}}\times \sum _{t=\mathrm{1}}^T{W}_t^{\mathrm{s},\mathrm{lw}},\text{(42)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{2}}^{\mathrm{lw}}=c_{\mathrm{2}}^{\mathrm{lw}}\times \sum _{t=\mathrm{1}}^T{W}_t^{\mathrm{s},\mathrm{lw}}\text{(43)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{3}}^{\mathrm{lw}}=c_{\mathrm{3}}^{\mathrm{lw}}\times \sum _{j=\mathrm{1}}^J\sum _{t=\mathrm{1}}^T{\displaystyle \frac{Q_{t,j}^{\mathrm{p},\mathrm{lw}}\times P_j^{\mathrm{lw}}}{Q_{j,\max }^{\mathrm{p},\mathrm{lw}}}},\end{array}$$

   where $M_{\mathrm{c},\mathrm{1}}^{\mathrm{lw}}$, $M_{\mathrm{c},\mathrm{2}}^{\mathrm{lw}}$, and $M_{\mathrm{c},\mathrm{3}}^{\mathrm{lw}}$ represent the water resource fees paid to the government, water fees paid to reservoir managers, and the electricity fees in Zhoushan Islands, respectively (RMB); $c_{\mathrm{1}}^{\mathrm{lw}}$, $c_{\mathrm{2}}^{\mathrm{lw}}$, and $c_{\mathrm{3}}^{\mathrm{lw}}$ denote the constant vectors, representing the unit price of water resources, water, and electricity in Zhoushan Islands, respectively (RMB per m³). Δt is the time step; i is the index of a reservoir; j is the index of a pumping station; I denotes the number of reservoirs; J denotes the number of pumping stations in Zhoushan Islands; $W_t^{\mathrm{s},\mathrm{lw}}$ denotes the amount of water supply for plants using local reservoir water at the tth time step (m³); $P_j^{\mathrm{lw}}$ denotes the supporting motor power of the jth pumping station (Kw); $Q_{t,j}^{\mathrm{p},\mathrm{lw}}$ denotes the flow through the jth pumping station at the tth time step (m³ s⁻¹); $Q_{j,\max }^{\mathrm{p},\mathrm{lw}}$ denotes the upper flow boundary of the jth pumping station in Zhoushan Islands (m³ s⁻¹).

2. operating costs for water supply using imported water ($M_{\mathrm{c}}^{\mathrm{iw}}$, RMB),

   $$\begin{array}{}\text{(44)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c}}^{\mathrm{iw}}=M_{\mathrm{c},\mathrm{1}}^{\mathrm{iw}}+M_{\mathrm{c},\mathrm{2}}^{\mathrm{iw}}+M_{\mathrm{c},\mathrm{3}}^{iw},\text{(45)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{1}}^{\mathrm{iw}}=c_{\mathrm{1}}^{\mathrm{iw}}\times \sum _{t=\mathrm{1}}^T{W}_t^{\mathrm{s},\mathrm{iw}},\end{array}$$
   $$\begin{array}{}\text{(46)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{2}}^{\mathrm{iw}}=c_{\mathrm{2}}^{\mathrm{iw}}\times \sum _{t=\mathrm{1}}^T{W}_t^{\mathrm{s},\mathrm{iw}},\text{(47)}& {\displaystyle }& {\displaystyle }{M}_{\mathrm{c},\mathrm{3}}^{\mathrm{iw}}=c_{\mathrm{3}}^{\mathrm{iw}}\times \sum _{j=\mathrm{1}}^J\sum _{t=\mathrm{1}}^T{\displaystyle \frac{L_j+Q_{t,j}^{\mathrm{p},\mathrm{iw}}}{Q_{j,\max }^{\mathrm{p},\mathrm{iw}}}},\end{array}$$

   where $M_{\mathrm{c},\mathrm{1}}^{\mathrm{iw}}$, $M_{\mathrm{c},\mathrm{2}}^{\mathrm{iw}}$, and $M_{\mathrm{c},\mathrm{3}}^{\mathrm{iw}}$ represent the water resources fees paid to the government, water fees paid to the river managers, and electricity fees in the city of Ningbo, respectively (RMB); $c_{\mathrm{1}}^{\mathrm{iw}}$, $c_{\mathrm{2}}^{\mathrm{iw}}$, and $c_{\mathrm{3}}^{\mathrm{iw}}$ denote the constant vectors, representing the unit price of water resources, water, and electricity in the city of Ningbo, respectively (RMB/m³); $W_t^{s,iw}$ denotes the amount of water supply for plants using imported water at the tth time step (m³); $Q_{t,j}^{\mathrm{p},\mathrm{iw}}$ denotes the flow through the jth pumping station at the tth time step, and J is the number of pumping stations transferring water from the city of Ningbo, J=2, $Q_{t,\mathrm{1}}^{\mathrm{p},\mathrm{iw}}=Q_{t,\mathrm{2}}^{\mathrm{p},\mathrm{iw}}$. L_j denotes the length of the continental diversion pipeline using the jth pumping station (m), and $Q_{i,\max }^{\mathrm{p},\mathrm{iw}}$ denotes the upper flow boundary of the ith pumping station for water transfer (m³ s⁻¹).

3. revenues (M_r, RMB),

   $$\begin{array}{}\text{(48)}& M_{\mathrm{r}}=b\times \left(\sum _{t=\mathrm{1}}^T{W}_t^{\mathrm{s},\mathrm{db}}+W_t^{\mathrm{s},\mathrm{th}}\right),\end{array}$$

   where b denotes the unit price of water supply revenue (RMB per m³).

The optimization model is subject to the following constraints:

$$\begin{array}{}\text{(49)}& {\displaystyle }& {\displaystyle }\left(\mathrm{1}\right)\mathrm{water}\mathrm{balance}:V_{t+\mathrm{1},i}=V_{t,i}+\left(I_{t,i}-Q_{t,i}^r\right)\mathrm{\Delta }t,\text{(50)}& {\displaystyle }& {\displaystyle }\left(\mathrm{2}\right)\mathrm{reservoir}\mathrm{storage}\mathrm{limits}:V_{t,i,\min }\le V_{t,i}\le V_{t,i,\max },\text{(51)}& {\displaystyle }& {\displaystyle }\left(\mathrm{3}\right)\mathrm{reservoir}\mathrm{release}\mathrm{limits}:Q_{t,i}^r\le Q_{t,i,\max }^r,\text{(52)}& {\displaystyle }& {\displaystyle }\left(\mathrm{4}\right)\mathrm{pumping}\mathrm{station}\mathrm{limits}:Q_{t,j}^p\le Q_{\max ,j}^p,\end{array}$$

where I_t,i is the inflow of the ith reservoir at the tth time step (m³ s⁻¹); V_t,i is the storage of ith reservoir at the tth time step (m³); V_min and V_max are the lower and upper storage boundaries, respectively (m³); $Q_{t,i,\max }^r$ is the maximum release of the ith reservoir at the tth time step (m³ s⁻¹). In some cases, $Q_{t,j}^p$ obtained by the RBF policies can be greater than $Q_{j,\max }^p$, and we will do the following step to modify $Q_{t,n}^r$.

$$\begin{array}{}\text{(53)}& {Q^{\prime }}_{t,n}^r={\displaystyle \frac{Q_{t,n}^r}{\sum _n^{N_j}{Q}_{t,n}^r}}\times Q_{j,\max }^p\end{array}$$

3.3 Model development

In this study, five input combination scenarios are considered to investigate whether the use of data-driven methods with climate forcing is efficient in inflow forecasts or not. These scenarios are described in Table 2. P_a represents antecedent precipitation, E_a represents antecedent evaporation, Q_a represents antecedent streamflow, P_f represents forecast precipitation, and E_f represents forecast evaporation.

Table 2 Five input combination scenarios.

Several strategies have been proposed in the literature to tackle a multi-step-ahead forecast task (Kline, 2004), such as the recursive, direct combination of direct and recursive strategies. In this study, we chose one of the most carried out strategies, i.e., the direct strategy (Ben Taieb et al., 2012), to forecast multi-step streamflow over the short-term horizon (1–7 d). In this case, the streamflow is forecasted using the following equations, using S3 as an example.

$$\begin{array}{}\text{(54)}& \begin{array}{rl}\mathrm{1}d:Q_{t+\mathrm{1}}^{\mathrm{f}}=& f\left(Q_t,Q_{t-\mathrm{1}},\mathrm{\dots },Q_{t-k},E_t,E_{t-\mathrm{1}},\mathrm{\dots },\right.\\ & \left.E_{t-k},P_{t-\mathrm{1}},\mathrm{\dots },P_{t-k}\right)\\ \mathrm{2}d:Q_{t+\mathrm{2}}^{\mathrm{f}}=& f\left(Q_t,Q_{t-\mathrm{1}},\mathrm{\dots },Q_{t-k},E_t,E_{t-\mathrm{1}},\mathrm{\dots },\right.\\ & \left.E_{t-k},P_{t-\mathrm{1}},\mathrm{\dots },P_{t-k}\right)\\ & \mathrm{\dots }\\ \mathrm{7}d:Q_{t+\mathrm{7}}^{\mathrm{f}}=& f\left(Q_t,Q_{t-\mathrm{1}},\mathrm{\dots },Q_{t-k},E_t,E_{t-\mathrm{1}},\mathrm{\dots },\right.\\ & \left.E_{t-k},P_{t-\mathrm{1}},\mathrm{\dots },P_{t-k}\right),\end{array}\end{array}$$

where f(⋅) is the mapping function between inputs and outputs, which can be modeled by LSTM, GRU, and GWO-LSSVM in our case. The hydrological variables normalized to the same scale of [0, 1] are used as the inputs in the three ML methods. The normalization equation is given as follows:

$$\begin{array}{}\text{(55)}& x^{\prime }={\displaystyle \frac{x-x_{\min }}{x_{\max }-x_{\min }}},\end{array}$$

where x and x^′ are the original and normalized values, respectively. x_min and x_max are the minimum and maximum values of the original series, respectively.

An issue with the ML methods is that they can easily overfit training data. To avoid this, the entire data are divided into three subsets in RNNs: (i) a training set, which is used to compute the gradient and update the weights and biases of the network; (ii) a validation set, over which the errors are monitored during the training process and is used to decide when to stop training; and (iii) a test set, which is used to assess the expected performance in the future. In addition, dropout is a regularization method where input and recurrent connections to LSTM and GRU units are probabilistically excluded from activation and weight updates while training a network. The strategies mentioned above have the effect of reducing overfitting and improving model performance in RNNs. Both LSTM and GRU are trained based on truncated backpropagation through time (BPTT) (Cheng et al., 2020), which uses a backpropagation network to update the parameters in iterations. The NSE function is used as the loss function to calibrate the LSTM and GRU models. As for LSSVM, we avoid overfitting by minimizing the NSE during the calibration and validation periods, while the test period is also used to assess the performance. In this study, January 2002 to December 2006 is used as the training period, while the validation and tests extend from January–December 2007 and January–December 2008, respectively.

The multi-reservoir operation optimization using inflow forecasts is performed over 1 year (1 April 2007–31 March 2008) with 25 reservoirs. The period is selected to ensure that it does not cover the calibration datasets. For the short-term forecasting and reservoir operation purpose, a forecast horizon of 1–7 d ahead is chosen. In this study, we use the parameterized MORDM approach to design operating policies for the multi-objective reservoir operations under uncertainty. The optimized operations are regulated based on both deterministic and uncertain forecast inflow. To keep it fair, we perform a simulation to generate deterministic and observed ensemble forecasts, which are each repeated 900 times. Using the uncertain streamflow forecasts (BMA, deterministic or observed ensemble forecasts) as policy inputs in the parameterized MORDM method, we can generate alternative RBF policies that are subject to all the constraints, and the objectives are evaluated over stochastic inflows. Under the parameterized MORDM, the decision variables in the optimization problem are not the volumes of water to be transferred from the city of Ningbo and the remaining 24 reservoirs each day. Instead, the decision variables are the parameters of the RBF policies. The best operation is obtained by conditioning the operating policies upon the following two input variables, e.g., the initial forebay water level and current inflow of reservoir. The optimization is solved at each time step (a particular forecast horizon, e.g., 1–7 d) by applying NSGA-II to search the space of decision variables and identify the islands' water allocation trajectories.

3.4 Results and discussion

3.4.1 Multi-step deterministic forecasts based on ML methods

We consider five different input scenarios described in Sect. 3.3. Table 3 demonstrates the forecast analysis carried out with the different configurations (input combination and forecast model), tabulating the NSE ranges for lead times from 1 to 7 d ahead over all reservoirs during the calibration, validation, and test periods. It can be seen that S1 using only the flow variables and S2 using only the antecedent climate variables are inferior to the other scenarios. The performance is generally improved when the flow variables are used in combination with the antecedent precipitation and evaporation under S3. However, in this case, the antecedent variables succeed in forecasting only 1 d ahead. The forecast performance decreases significantly as the forecast horizon increases from 1 to 7 d ahead. Herein, we suppose that the following precipitation and evaporation have been forecasted. It is clear that S4 and S5, with the forecast climate variables, make significant increments in streamflow forecasting. The NSE can remain relatively stable at different horizons. There are no apparent differences between the three forecast models during the calibration period. However, the two RNNs perform better than GWO-LSSVM during the validation period, while GWO-LSSVM outperforms during the test periods. Besides, given that GRU has more superficial structures and fewer parameters and requires less time for model training, it may be the preferred method for short-term streamflow forecast compared with LSTM. The same results have been obtained in Gao et al. (2020) when they used LSTM and GRU to model short-term rainfall–runoff relationships.

Table 3 NSE ranges ([min, max]) for all reservoirs with the different configurations during the calibration, validation, and test periods.

We aim to compare how the forecasted climate variables impact the streamflow forecast and reservoir operation performance. For the sake of brevity, S3 and S5 are compared in detail in the following section. Recall that S3 uses flow variables, antecedent precipitation, and evaporation as inputs, while S5 uses flow variables as well as the antecedent and forecast climate forcing. After assessing model validity, the next step is to compare the performance across the 24 reservoirs. The coefficient of variation (COV), defined as the ratio of the standard deviation of the inflow time series, is used to capture the varying characteristics of the incoming flow into the reservoir. Figure 6 reveals a strong negative relationship between COV and forecast performance under S3 at all lead times. The forecast performance decreases as the COV increases for all forecast models. This indicates that the more variation the flow has, the harder it is for data-driven methods to learn the flow pattern when there is not enough input information. However, the negative signal under S5 (Fig. 7) with forecasted climate variables (precipitation and evaporation in this study) is not as strong as that under S3, indicating that the forecast climate variables can help AI-based models' mapping functions between inputs and outputs. The improvements are more significant for the two RNN models, i.e., LSTM and GRU, than LSSVM. This result demonstrates that the efficiency of deep-learning RNN methods is better and more accurate than LSSVM.

Figure 6 NSE values at lead times of 1 to 7 d plotted against the coefficient of variation (COV) for all the 24 reservoirs during the period of (a) calibration, (b) validation, and (c) test under S3.

Figure 7 NSE values at lead times of 1 to 7 d plotted against the coefficient of variation (COV) for all the 24 reservoirs during the period of (a) calibration, (b) validation, and (c) test S5.

3.4.2 Multi-step stochastic forecasts based on BMA method

Based on the forecast results of three data-driven models in the calibration period, the BMA method determines weights for LSTM, GRU, and GWO-LSSVM models. The weights reflecting the performance of the ensemble models during the calibration period are shown only for lead times of 1 and 7 d for the sake of brevity under S3 and S5 in Fig. 8. The model weights reflect the comparative importance of all the competitive modeling predictions on one level. Figure 8 indicates that it is difficult to conclude which individual model provides the best prediction. For example, GRU outperforms the remaining two models for Hongqiao Reservoir, while LSTM performs best for Cenggang Reservoir in Fig. 8a1. Similar results can be obtained from Fig. 8b1. Comparatively, Fig. 8a2 shows that LSTM and GWO-LSSVM influence the BMA model more than GRU. This higher weight is assigned because the forecasts are more similar to observations than those less similar to observations using the BMA posterior processor. However, observed from Fig. 8b2, the prediction accuracy of GWO-LSSVM is seriously affected and much less than that of GRU. It is consistent with the results obtained in Fig. 7, indicating that RNNs outperform GWO-LSSVM when there is more input information under S5. Overall, model uncertainty always exists whether forecast climate variables are involved or not, and it is necessary to analyze and evaluate the model uncertainty involved using BMA.

Figure 8 Weights of three individual forecast models for the BMA model for all reservoirs under (a) S3 and (b) S5.

To access model validity, the evaluation of the modeled streamflow is performed over calibration, validation, and test periods using NSE, RMSE, and MAE metrics. Table 4 shows the performance metric ranges for all 24 reservoirs of BMA methods under S3 and S5. Apparently, both the replicative (forecast performance in calibration sets) and predictive (forecast performance in validation and test sets) validity under S5 for forecast horizons are significantly better than those under S3. For example, Fig. 9 demonstrates the improvement rates in terms of NSE, RMSE, and MAE of the BMA model compared with the three individual models. BMA produces the maximum NSE, minimum RMSE, and minimum MAE during the calibration period for both two scenarios, indicating that BMA has the best goodness of fit. This is because the weights are derived according to the individual forecast model in this period. With respect to validation and test periods, the BMA method shows better forecasts than the three comparative models except for the GRU modeling validation datasets under S5. Thus, it is shown that the BMA model matches the actual streamflow well.

Figure 9 Improvement rates in terms of averaged (a) NSE, (b) RMSE, and (c) MAE of the BMA model for forecasts as compared with the three individual models.

Table 4 Performance metric ranges ([min, max]) for all 24 reservoirs of BMA methods under S3 and S5.

The model validity is then assessed using (i) hydrographs and (ii) scatter plots of observed and modeled streamflow, as shown in Figs. 10 and 11. Herein, we only show three reservoirs, i.e., Hongqiao (the largest reservoir), Goushan (the medium reservoir), and Nanao (the smallest reservoir), for the sake of brevity. From Fig. 10, it is clearly shown that the modeled streamflow deviates gradually from the 1:1 line, and the forecast skill decreases with the increase of lead time under S3 as expected, which is consistent with the statistical results shown in Table 4. In contrast, the scatters of the observed and modeled streamflow implemented with forecasted climate variables fit well across the 1:1 line at different lead times under S5, observed from Fig. 11. The performance for Hongqiao Reservoir is affected explicitly by an extreme peak event that hit the reservoir during the calibration period shown in Fig. 10, which does not occur over the training set of data. This causes heavy underestimations in the streamflow forecast. A more extended calibration period is required to improve the performance over such extreme peak flow events. However, the BMA method performs well on this extreme peak flow in Hongqiao Reservoir at all lead times when the forecast climate forcing is applied as inputs. This is because the reservoirs in Zhoushan Islands have relatively small drainage areas, and thus the flow is concentrated in a very short time after an extreme rain event.

Figure 10 Forecast results of (a) Hongqiao, (b) Goushan, and (c) Nanao reservoirs under S3.

Figure 11 Forecast results of (a) Hongqiao, (b) Goushan, and (c) Nanao reservoirs under S5.

We use the Monte Carlo simulation method to generate BMA ensemble forecasts. The number of simulations is set as 1000 in this study. To demonstrate the optimization results of multi-reservoir operations based on the data-driven forecast models under uncertainty, 90 % confidence intervals associated with the deterministic predictions at BMA are further calculated. The confidence interval provides more alternatives that are possibly useful for a tradeoff between multiple objectives, such as flood control, hydropower generation, and improved navigation (Zhang et al., 2015). The interval performance metrics of Cr and D described in Sect. 2.3 are adopted to assess the performance of uncertain forecasts. Table 5 displays the averaged metrics for all the 24 reservoirs under S3 and S5. Both indicators under S5 are superior to those under S3. The 90 % streamflow interval between the 5th and 95th percentiles of some representative reservoirs, e.g., Hongqiao, Goushan, and Nanao reservoirs, are presented in Figs. 12 and 13. The results are consistent with those in Figs. 10 and 11. It is observed from Fig. 12 that the streamflow interval fails to capture the extreme peak flow for Hongqiao Reservoir under S3. The BMA performs gradually worse with increasing lead times for the three reservoirs. However, in Fig. 13, the red dots represent the observed streamflow, most of which are covered by the 90 % interval at both 1 and 7 d ahead. Therefore, the forecast climate variables will be conducive to reducing the predictive uncertainty of real-time streamflow forecasting.

Figure 12 The 90 % streamflow interval of the BMA method under S3.

Figure 13 The 90 % streamflow interval of the BMA method under S5.

Table 5 Ranges of interval performance metrics ([min, max]) for all the 24 reservoirs under S3 and S5.

3.4.3 Multi-objective reservoir operation performance evaluation

The optimized operations are regulated based on both deterministic and uncertain forecast inflow. To demonstrate the relationship between the conflicting objectives, a set of Pareto solutions over a 7 d horizon at different periods under S5 is given as an example, as shown in Fig. 14. The optimization using the Pareto concept allows the operator to choose an appropriate solution depending on the prevailing circumstances and analyze the tradeoff between the conflicting objectives. In each of the plots, the water deficiency ratio of Daobei Plant and the sum of the remaining plants are plotted on the x and y axes, respectively. The color of the markers indicates the net operating costs, with colors ranging from red, representing low value, to blue, representing high value. Thus, an ideal solution should be located at the left corner (low value of the water deficiency ratio of Daobei Plant and the sum of the remaining three plants) of the plot and represented by a red (low net operating costs) marker. The black arrows have been added in the figure to guide the reader in understanding the directions of optimization. Generally, the water deficiency ratio of Daobei Plant has an inverse relationship with that of the sum of the remaining plants (inverse relationship; i.e., the former decreases with the increase of the latter). In contrast, the water deficiency ratio of the remaining three plants has a positive relationship with the net costs (positive relationship; i.e., the former increases with the increase of the latter).

Figure 14 A set of Pareto solutions at different periods over a 7 d horizon under (a) deterministic and (b) uncertain forecasts.

It is interesting to compare the performances associated with deterministic and uncertain forecasts. Uncertain conditions (Fig. 14b) show a much broader scale on the three objectives than deterministic conditions (Fig. 14a). For instance, uncertain forecasts produce the water deficiency ratio of Daobei Plant, ranging from −40 % to 80 %, during 12 to 18 August 2007, while deterministic forecasts have a smaller range, with a value from 30 % to 100 %. The water supply deficits under deterministic forecasts are due to the high demand happening in August, which can be mitigated when informing the reservoir operations with uncertain forecasts. In this way, we expect that if the ensemble streamflow forecasts are used in a stochastic optimization scheme, the reservoir operation could be further enhanced because the optimization considers possible uncertainty provided by uncertain forecasts and thus takes advantage of correcting the influences of uncertainty.

1. Performance evaluation with different forecast skills. In general, forecasts are always useful for reservoir operations. The annual revenues, costs, and water supply reliability are chosen as metrics to compare the performance of the operating policies derived from different configurations. Reliability is a measure of how well the water demand for users is met in a water transfer system. In this case, reliability is expressed as a percentage. The system performances are averaged over a set of solutions. The annual values during the period from 1 April 2007 to 31 March 2008 at various configurations are provided in Table 6 with two decision horizons of 1 and 7 d. The multi-reservoir operation based on observation is designed as a benchmark. It can be seen from Table 6 that the performance indicators from the 1 d forecast horizon are better than those from 7 d using deterministic inflows (in the case of observed and forecasted inflows). Two scenarios (S3 and S5) with the 1 d forecast horizon show similar operating performance, which is consistent with the performance of the inflow forecast listed in Table 3. Recall again that S3 uses flow variables, antecedent precipitation, and evaporation as forecast inputs, while S5 uses flow variables as well as the antecedent and forecast climate forcing. In contrast to S3, the operating results of S5 with a 7 d forecast horizon are closest to that of the observation. This is due to the improved inflow forecast performance under S5. However, it is depicted in Table 6 that the indicator of water supply reliability and net costs under S5 are inferior to those under S3. As for the stochastic forecasts, S5 outperforms S3 with lower net costs and approximate water supply reliability. In this case, the improved performance may not lead to improved decisions in deterministic forecasts.

   The results obtained in Table 6 show that system performance derived from the observed inflows is inferior to that from other configurations. This finding cannot confirm the effectiveness of inflow forecasts. The reason for that is the forecast inflows may overestimate the actual inflows. For example, the mean value (0.14 m³ s⁻¹) of the observed inflow of Hongqiao Reservoir is lower than that of the forecasted inflow (0.17 m³ s⁻¹). In this case, the good performance presented in Table 6 is “fake”. That is to say, although decision makers can follow the strategies determined by the forecasted inflows, the system performance should be assessed using the actual inflows (i.e., observed inflows). We further re-evaluate the operating strategies optimized from different configurations mentioned above using the observed inflows. The performance metrics are listed in Table 7. It is expected that the results can reveal the maximum efficiency and reliability that could be achieved based on accurate information. In general, the indicator values under deterministic forecasts in Table 7 are reduced compared with those in Table 6. The reason is that reservoir operating decisions in Table 6 are optimized based on a higher inflow series.

   In terms of both deterministic and uncertain forecasts, net operating costs of S5 are improved significantly compared with that of S3, while water supply reliability is increased slightly. This result may suggest that improved forecasts are more skillful in making decisions when using forecast climate variables as inputs. We highlight that this result we obtained is specific to the Zhoushan Islands. Indeed, many studies show that higher forecast performance did not lead to better operation decisions (Chiew et al., 2003; Goddard et al., 2010; Turner et al., 2017). However, some researchers draw the same conclusions as us. For instance, Anghileri et al. (2016) declared that inflow forecasts with accurate weather components would produce much smaller water supply deficits. Moreover, Anghileri et al. (2019) found that preprocessed forecasts (higher performance) were more valuable than the raw forecasts (less performance) regarding two operation performance metrics, i.e., mean annual revenues and spilled water volume.

   There is also an interesting finding that the operating performance upon deterministic forecasts deteriorates, while the performance upon uncertain forecasts can stay relatively stable. This implies that the use of uncertain forecasts in reservoir operation can be more efficient and reliable than that of deterministic forecasts. The reason is that in a stochastic optimization scheme, the value can be further enhanced because the optimization can account for the total uncertainty provided by the ensemble forecasts. Similar results were obtained by Roulston and Smith (2003), who reported that the hydroelectric power production derived from the ensemble forecasts was increased compared with the deterministic forecasts. Boucher et al. (2012) also found that stochastic forecasts outperformed deterministic ones with the lower turbinate flow, higher generation production, and less spillage during a flood period. Overall, in most cases, a noticeable improvement can be achieved through the use of the stochastic decision-making assistance tool.

   We then assess the performance metrics of water supply reliability over different seasons. It is noted in Fig. 15 that the deterministic forecasts are less skillful than the uncertain forecast when used in spring (JFM), summer (AMJ), autumn (JAS), and winter (OND) with the two forecast horizons. Although the operating performance using the deterministic forecast is lower due to its deterministic character, the main characteristics of the relationship between the forecast quality and value remain unchanged. That is to say, the benefits of considering the forecasts are more significant when the forecast quality is higher. It indicates that the optimization is capable of exploiting efficient information to improve reservoir operations. In our multi-objective optimization modeling, we would like to make the best use of water resources and maximize water supply. However, the operating performance in autumn shows a lower value with respect to that in other seasons. This is because the water demand in autumn is usually much higher. The shortage does not imply the non-effectiveness of our proposed forecast-based management framework but is due to the limitation of available water and pumping station capacity.

2. Performance evaluation with different forecast horizons. The impact of different forecast horizons on the operation performance is further evaluated under different configurations, as shown in Fig. 16. It is noted that the operating policy optimized from uncertain forecast inflows upon S5 outperforms that from S3. In terms of deterministic conditions, S5 improves the operation on the metrics of water supply reliability of Daobei Plant, water supply reliability of the other plants, and net costs with a variation of 2.11 %–13.58 %, 2.74 %–7.38 %, and −19.94 % to −10.30 %, respectively, compared with S3. As for uncertain conditions, S5 improves by 0.24 %–1.90 %, 0.06 %–1.32 %, and −59.45 % to −176.19 %, respectively. Although the increments in water supply reliability are not insignificant, S5 can secure water demand with much lower operating costs than S3, which decision makers value most. Furthermore, uncertain forecasts produce an improved ratio of 31.52 %–65.01 %, 19.98 %–46.60 %, and −116.45 % to −56.95 % than deterministic forecasts regarding the three metrics, respectively. Our results again highlight that uncertain forecasts are more valuable than deterministic forecasts when designing forecast-informed reservoir operations.

   With an increase in forecast horizon from 1 to 7 d, the performance in water supply reliability and net operating costs upon deterministic conditions are generally reduced. This suggests that considering a longer forecast horizon (up to 7 d) does not necessarily improve reservoir operation without future forecast climate variables as inputs (low forecast quality). The reduced performance in water supply reliability might be due to the fact that the optimization explores strategies to secure the whole water demand in a longer horizon, which results in a sacrifice in reliability on some particular days. This result is similar to the finding proposed in Xu et al. (2014), who argued that the use of longer horizon (an efficient forecast horizon longer than 1 d) inflows could not improve hydropower performance when they set the forecast horizon as 1–5 d. Nevertheless, the increasing forecast horizon may not generate improved or decreased water supply reliability in uncertain conditions. Approximate water supply volume can lead to similar revenues or fees paid to the government and managers (water fees and water resources fees). Accordingly, the growing trend in net costs is caused by the increased operating costs, mainly dominated by electricity prices, when the multi-reservoir is operated to supply the demand in a longer horizon. In this case, the operation performance varies at different conditions. This demonstrates that the relationship between the forecast horizon and reservoir operation is rather complex and depends not only on the configurations (i.e., inflow forecast quality and uncertainty) used to determine operating rules, but also on the performance metrics used to assess operation.

Figure 15 Seasonal system performance of water supply reliability.

Figure 16 Annual system performance with different forecast horizons.

Table 6 Annual system performance using forecast inflow information.

Table 7 Annual system performance using observed inflow information.

4 Limitations and future work

Our work suffers from some limitations which could be overcome in future studies. One of the limitations is that only one single indicator was used to calibrate the forecast models, while multiple indicators were used in assessing the performance of the models. It would be a more fair practice to use multi-criteria to do both calibration and assessment, and this could be interesting for future work. Another limitation is that we used the average observed price to calculate the revenues and operating costs. In an operational and deregulated market setting, the prices may fluctuate significantly (Anghileri et al., 2019). For instance, forecasting electricity prices is likely to improve short-term operation efficiency significantly. The combined effects of price and streamflow forecasts on water resource allocation are worth investigating in future studies. Our study also suffers from the drawback that instead of using the short-term weather forecasts from the Global Forecast System (GFS) or European Centre for Medium-Range Weather Forecasts (ECMWF) model (Choong and El-Shafie, 2015; Schwanenberg et al., 2015; Peng et al., 2018; Ahmad and Hossain, 2019; Liu et al., 2019), we used the observed weather conditions as alternatives, which may result in an overestimation in forecast quality. However, forecast uncertainty and error generally grow with lead time. The usefulness of the forecast information is reduced with the increase of the forecast horizon and thus the operating performance. This may influence the finding we highlight above that the relationship between the forecast horizon and reservoir operation is not constant and specific. It would be interesting to analyze the reservoir operation performance when accounting for an ensemble numerical weather prediction.

5 Conclusions

In this study, we proposed an AI-based management methodology to assess forecast quality and forecast-informed reservoir operation performance together. The approach was tested on a water resources allocation system in Zhoushan Islands, China. Specifically, the findings are summarized below.

A data-driven reservoir inflow forecasting system using ML methods (LSTM, GRU, and GWO-LSSVM) was first developed with a comprehensive calibration–validation–testing framework. The validity of the deterministic forecast was demonstrated by applying it over 25 reservoirs with varying climate and hydrological characteristics. Results showed that the more variation the streamflow has (a high COV value), the harder it was for the ML methods to learn the flow pattern when there was not enough input information. The forecast skill deteriorated with increasing lead times under such scenarios. However, short-term forecast climate forcing was efficient and scalable in forecasting the multi-reservoir inflow over the forecast horizon (1–7 d). LSTM and GRU models generated comparable performance under different configurations. Given that GRU has simpler structures and fewer parameters and required less time for modeling, it might be the preferred method for streamflow forecasts than LSTM.

Then we used BMA to generate stochastic inflow scenarios for quantifying uncertainty based on LSTM, GRU, and GWO-LSSVM deterministic forecasts. The results demonstrated that it was difficult to conclude which individual model provided the best prediction, but the BMA did display better forecast skills in comparison to the individual ones. Including one scenario with antecedent conditions and one scenario with both antecedent and forecast information, two input combination scenarios were compared on the uncertain forecast performance in detail. The comparison indicated that forecast climate variables would help reduce the predictive uncertainty of short-term streamflow forecasting.

The forecasting scheme was further coupled with a multi-objective reservoir operation model to optimize water resources allocation. Using a MORDM approach, we identified strategies that were useful for a tradeoff between water supply reliability and operating costs in Zhoushan Islands. A rolling horizon scheme was employed to obtain an optimal operating policy over the horizon of 1–7 d. The long-term assessment over a year based on deterministic and stochastic forecasts showed quite different performances in terms of water supply reliability and net operating costs. Our averaged annual results showed that uncertain forecasts were more valuable than deterministic forecasts. The operating benefits of considering the forecasts were more significant when the forecast quality was higher. Similar results could be obtained at a seasonal scale. While showing the unquestionable benefit of implementing forecast-based reservoir operations, our results also demonstrated that the relationship between the forecast horizon and reservoir operation was complex and depended on the operating configurations (forecast quality and uncertainty) and performance measures for the Zhoushan Islands system.

Overall, the developed AI-based management framework has demonstrated a clear advantage in quantifying the uncertainties of inflow forecasts to improve the overall system performance of water allocation systems. Such a framework can be further applied to other study sites with similar problems. However, the results we obtained in this study are only specific to the Zhoushan Islands and should be applied to other study sites with care.