Combining Grey Relational Analysis and a Bayesian Model Averaging Method to Derive Monthly Optimal Operating Rules for a Hydropower Reservoir

Various regression models are currently applied to derive functional forms of operating rules for hydropower reservoirs. It is necessary to analyze and evaluate the model selecting uncertainty involved in reservoir operating rules for efficient hydropower generation. Moreover, selecting the optimal input variables from a large number of candidates to characterize an output variable can lead to a more accurate operation simulation. Therefore, this paper combined the Grey Relational Analysis (GRA) method and the Bayesian Model Averaging (BMA) method to select input variables and derive the monthly optimal operating rules for a hydropower reservoir. The monthly input variables were first filtered according to the relationship between the preselected output and input variables based on the reservoir optimal deterministic trajectory using GRA. Three models, Particle Swarm Optimization-Least Squares Support Vector Machine (PSO-LSSVM), Adaptive Neural Fuzzy Inference System (ANFIS), and Multiple Linear Regression Analysis (MLRA) model, were further implemented to derive individual monthly operating rules. BMA was applied to determine the final monthly operating rules by analyzing the uncertainty of selecting individual models with different weights. A case study of Xinanjiang Reservoir in China shows that the combination of the two methods can achieve high-efficiency hydropower generation and optimal utilization of water resources.


Introduction
Hydropower is a clean and renewable energy source and accounts for 20% of electricity generation worldwide [1].Considering the economic, technical, and environmental benefits of hydropower, most countries, especially for developing countries, usually have a tremendous and ever-intensifying need for electricity, and they also possess the most significant remaining hydropower potential [2].A hydropower reservoir is one of the efficient ways to explore the electricity generation reliability and economic benefit [3].Accordingly, developing hydropower generation by implementing feasible reservoir operation rules is being recognized as a strategic issue [4].The conventional hydropower reservoir operations currently prescribe reservoir releases based on limited criteria such as current storage level and inflow [5].Compared with the conventional operation, the deterministic optimization model [6] of hydropower generation can detect an optimal solution for better utilization of available resources with the increasing complexity and interdependency of systems in reservoir management [7].However, the deterministic optimization with perfect inflows or other system inputs is challenging to Water 2018, 10, 1099 2 of 20 apply to real operations.Reservoir operation decisions should be made as needed in "real-time" with the limited foresight of future conditions [8].
There are two commonly considered forms of reservoir operating rules, namely the reservoir operating function and the operating chart [9].The reservoir operating function directly derives the operating rules from a deterministic optimization model that can better inherit the optimization efficiency in hydropower generation than the operating chart.Various functional forms are successfully applied to the deduction of reservoir operating rules, including Multiple Linear Regression Analysis (MLRA) methods [10], artificial neural network (ANN) approaches [11], adaptive-network-based inference system (ANFIS) techniques [12], support vector machine (SVM) method [13], etc.Although the above operating rules have been successfully used, except the MLRA model, the other machine learning models are black box models, and all models will have model bias.Therefore, it is necessary to analyze and evaluate the model selecting uncertainty involved in reservoir operating rules.The Bayesian model averaging (BMA) method can overcome this uncertainty of selecting models by conditioning, not on a single "best" model, but on the entire ensemble of statistical models first considered [14,15].The BMA method is widely applied in many fields, such as hydrological forecasting [15], social and health sciences [16], stroke risk [17], and groundwater modeling [18,19].However, rare applications of this technique in the area of deduction of reservoir operating rules have been recently reported [20].
The reservoir operating rules determining the operation process at which the output variables, such as the final water level, release, or power output, are performed are based on the currently available reservoir information referred to as input variables.For a more accurate, parsimonious, and physically interpretable simulation for the reservoir operation process, selecting the optimal input factors from a large number of candidates to characterize an output variable is indispensable [21].Input variables selection approaches are commonly divided into two categories, model-based and model-free methods [22].Model-based methods can efficiently figure out the best inputs from different sets of inputs according to a lot of calibration and validation processes, but which are considered computationally intensive and time-consuming [23].Different from the model-based method, model-free methods work based on the relationship between the input and output variables, as measured by interclass distance, statistical dependence, or an information-theoretic measure [24,25].Nevertheless, few studies have concentrated on the selection of input factors using model-free methods while deriving the reservoir operating rules.Ji et al. [26] firstly implemented the multiple linear regression algorithms for selecting appropriate time, space, and energy factors to optimize the hydropower generation operating function for cascaded reservoirs.They [27] then applied the rough sets theory to remove redundant attributes of the input variable sets deriving the flood operation rules.Yang et al. [28] proposed a cascade-reservoir input-variables selection method considering the relations between input variables and decision-making in optimal reservoir operation using the extra-tree model.In addition, the relationship between the preselected output and input variables is not constant during the different periods.Thus, the monthly optimal input variables are needed for efficiently modelling of operating rules.
Grey relational analysis (GRA) proposed by Deng [29] is an impacting measurement method in grey system theory that analyzes uncertain relations between one main factor and all the other factors in a given system.It could be used to measure the approximate correlation between sequences with convenient and small procedures [30].The GRA method is successfully applied to select the optimal input variables in many fields [31][32][33], including agriculture, traffic, industrial engineering, education, but rare studies focus on reservoir operation problems [34,35].In this paper, GRA is first applied to derive the optimal input factors based on the relationship between the preselected output and output variables in the area of hydropower reservoir operation.
For efficient hydropower generation, the primary purpose of this study is to combine the GRA and BMA method to derive the monthly optimal operating rules for a hydropower reservoir.An optimal deterministic operation model of reservoir hydropower generation is first established and solved.
Water 2018, 10, 1099 3 of 20 Then, the optimal monthly input factors are determined based on the correlation between preselected output and input variables according to the optimal trajectory using GRA.PSO-LSSVM, ANFIS, and MLRA are further applied to derive individual operating rules.Lastly, BMA is used to determine the final reservoir operating rules.The flowchart of the study is shown in Figure 1.A case study of Xinanjiang Reservoir in China is applied to demonstrate the combination of the two methods.
optimal deterministic operation model of reservoir hydropower generation is first established and solved.Then, the optimal monthly input factors are determined based on the correlation between preselected output and input variables according to the optimal trajectory using GRA.PSO-LSSVM, ANFIS, and MLRA are further applied to derive individual operating rules.Lastly, BMA is used to determine the final reservoir operating rules.The flowchart of the study is shown in Figure 1.A case study of Xinanjiang Reservoir in China is applied to demonstrate the combination of the two methods.

Model Formulation
As this study aims to derive the optimal operating rules for a hydropower reservoir, two objectives while satisfying all kinds of constraints are considered here: one is to maximize the hydropower generation [36], and the other is to maximize the firm output.A maximin model [37] is selected to optimize the firm power, as shown in Equation ( 2), which has been proved to be capable of obtaining the same minimal output as the equal output method with more direct and reliable results.
( ) The equivalent transformation [37] is to maximize the hydropower generation with a penalty function for a violation of firm power output (Equation ( 3)).The reservoir optimization model operates by determining an optimal release or water level for a reservoir over the whole operation period.In this study, the monthly reservoir water level is set as the decision variable.The release of the reservoir is defined with a predefined water level according to the water balance equation: where J is the objective (KW•h); A is the power output coefficient of the hydropower station; Ht is the water head of the hydropower station, Ht = Zt -Zd,t, Zt is the reservoir water level, Zd,t is the tailwater level, (m); qt is the turbine release of hydropower station at time step t (m 3 /s); Nt is the power output of hydropower station at time step t (MW); Nf is the firm power of the hydropower reservoir (MW); Δt is the time step (s); and μ is the penalty factor, which is generally regarded as a large number.
The operation model of a reservoir is subject to the following constraints: a. Water storage constraint:

Model Formulation
As this study aims to derive the optimal operating rules for a hydropower reservoir, two objectives while satisfying all kinds of constraints are considered here: one is to maximize the hydropower generation [36], and the other is to maximize the firm output.A maximin model [37] is selected to optimize the firm power, as shown in Equation ( 2), which has been proved to be capable of obtaining the same minimal output as the equal output method with more direct and reliable results.
The equivalent transformation [37] is to maximize the hydropower generation with a penalty function for a violation of firm power output (Equation (3)).The reservoir optimization model operates by determining an optimal release or water level for a reservoir over the whole operation period.In this study, the monthly reservoir water level is set as the decision variable.The release of the reservoir is defined with a predefined water level according to the water balance equation: where J is the objective (KW•h); A is the power output coefficient of the hydropower station; H t is the water head of the hydropower station, H t = Z t -Z d,t , Z t is the reservoir water level, Z d,t is the tailwater level, (m); q t is the turbine release of hydropower station at time step t (m 3 /s); N t is the power output of hydropower station at time step t (MW); N f is the firm power of the hydropower reservoir (MW); ∆t is the time step (s); and µ is the penalty factor, which is generally regarded as a large number.
The operation model of a reservoir is subject to the following constraints: a.
Water storage constraint: Water 2018, 10, 1099 4 of 20 c. Turbine release constraint q min,t ≤ q t ≤ q max,t d.
Power output constraint N min,t ≤ N t ≤ N max,t (7) where t is the index of time step; Q t and U t are the inflow and outflow of reservoir at time step t (m 3 /s), respectively; S t is the storage of reservoir at time step t (m 3 ); and min and max are the lower and upper boundaries, respectively.In this study, a Dynamics Programming (DP) method was applied to solve the deterministic operation model to find the optimal or near-optimal trajectory (the monthly power output and reservoir water level).

Description of the Input and Output Variables
The power output N is selected as the output variable here as this paper aims to generate rules for hydropower generation.Input variables are distributed into two types, namely time and energy factors.The time factors representing the state of a reservoir includes the initial water level Z 0 , inflow Q, and maximum water level Z a , as in Equation ( 7) [27]: where φ is the relation function between the storage and the water level; V 0 is the initial reservoir storage at time step t (m); and T is the time step, T = 2.63 × 106 s.In this paper, the incoming hydropower E f (kW•h), the storing hydropower E s (kW•h), and the interaction between incoming and storing hydropower E fs (kW 2 •h 2 ) are considered as the energy factors, represented in Equation ( 9) [27].
where A indicates the output power coefficient a hydropower station; ϕ is the mapping function between the outflow and tailwater level; V t is the final reservoir storage when all inflows are impounded at time t step (m 3 ); and Z 0 represents the downstream water level when there is no outflow from a reservoir (m).

Grey Relational Analysis (GRA) Selecting the Monthly Impact Factors
The GRA method can investigate the uncertain relationships between one main factor and all the other factors in a system with five layers [38].
Step 1: Establish the grey relational matrix.
Step 2: Normalize the grey relational matrix.
A linear normalization of the grey relational matrix using Equation ( 12) is performed in the range from zero to one.
where X j = {X 0 ;X i }, X 0 is the output, X i is the i-th input; GRN j = {GRN 0 ; GRN i }, GRN j,t is the i-th grey relational normalization, GRN i is the grey relational normalization of the i-th input factor, GRN 0 is the grey relational normalization of output variables; max t X j and min t X j are the maximum and minimum values of X j ; Step 3: Calculate the grey relational coefficient.
where η i,t is the grey relational coefficient between GRN i,t and GRN 0,t ; ρ = 0.5 is taken for proper stability of outcomes with moderate effects.
Step 4: Calculate the grey relational grade.
where r i is the grey relational grade between X(i) and X( 0 ).
The grey relational rank is based on the grey relational grade.The experimental parameters of the higher ranked input variable are closer to the output variable.

Individual Regression Models
For individual models, consider a given training set of T data points {x i , f i } T i=1 with input data x and output f.

Particle Swarm Optimization-Least Squares Support Vector Machine Model
The LSSVM model is a non-linear regression forecasting method proposed by Suykens and Vandewalle [39].In future space, the SVM model takes the form f (x) = ω T ϕ(x) + b, where the nonlinear mapping ϕ( ) maps the input data in to a higher dimensional feature space, and b is the bias.Note that the dimension of ω is not specified.In LSSVM, for the function estimation, the following optimization problem is formulated: subject to the equality constraints: This corresponds to a form of ride regression.The Lagrangian is given b by y with Lagrange multipliers α k .
Water 2018, 10, 1099 This finally results in the following LSSVM model for the function estimation: where Ψ(x, x i ) is the kernel function.There are four commonly used kernel functions available: linear kernel, polynomial kernel, radial basis kernel function (RBF), and the sigmoid function [40].The RBF is widely used as it not only computationally simpler than other functions, it also maps the training data to an infinite-dimensional space in a nonlinear manner [41][42][43].Thus, this study selects the RBF as the kernel function.The RBF has two parameters, C and σ, which are calibrated using the particle swarm optimization (PSO) algorithm.The PSO is an evolutionary algorithm inspired by the feeding behavior characteristic of a bird flock [44,45].

Adaptive Neural Fuzzy Inference System Model
ANFIS is a combination of ANN and a fuzzy inference system (FIS).To obtain a better modelling system, ANN can be combined with FIS to improve speed, fault tolerance, and adaptiveness [46].The network structure is capable of adjusting the shape of the membership functions and the consequence parameters that form the fuzzy rules by minimizing the difference between the output and provided targets [47].To illustrate those two procedures of an ANFIS, for simplicity, two inputs, x 1 and x 2 , and one output, y, are assumed.For a first-order Sugeno fuzzy model, a typical rule set with two fuzzy if-then rules can be expressed as: where A and B are the fuzzy sets, and p i ,q i , and r i are linear parameters in the "then" part (consequent part) of the first-order Sugeno fuzzy model.ANFIS is a feed-forward neural network with five layers (Figure 2), and a brief introduction of the model is as follows: Layer 1: Every node in this layer is an adaptive node with a node output defined as: where x 1 and x 2 are the input nodes, A and B are the linguistic labels, and µ(x 1 ) and µ(x 2 ) are the membership functions, which are usually adopted a bell shape with maximum equal to 1 and minimum equal to 0, as follows: Thus, , where α i = {α i ,b i ,c i } are the premise parameters, and i is the index of linguistic label of each input variable.
Level 2: Every node in this layer is a fixed node label Π with the node function to be multiplied by input vectors to serve as an output.The output ω represents the firing strength of a rule.For instance: Level 3: Every node in this layer is a fixed node, marked by a circle and label N, with the node function to normalize the firing strength by calculating the ratio of the ith node firing strength to the sum of all rules' firing strength.
Water 2018, 10, 1099 7 of 20 Level 4: The fourth layer is called the implication layer.The consequence of each rule is calculated as a linear combination of the input variables, as described by Takagi and Sugeno [48], and then multiplied by its associated normalized firing strength: where i is the ith node's output from the previous layer, and {pi, qi, ri} are the consequence parameters that increase with the number of input variables.Level 5: In the fifth layer, all the incoming signals are summed to compute the simulated power output: For the data training, the hybrid learning algorithm of the ANFIS combines the gradient method with the least squares method to update the parameters [49].In the forward pass of the learning algorithm, consequent parameters are identified by the least squares estimate.In the backward pass, the error signals, which are the derivatives of the squared error with respect to each node output, propagate backward from the output layer to the input layer.In this backward pass, the premise parameters are updated by the gradient descent algorithm.
Water 2018, 10, x 7 of 20 where i ϖ is the ith node's output from the previous layer, and { , , } i i i p q r are the consequence parameters that increase with the number of input variables.Level 5: In the fifth layer, all the incoming signals are summed to compute the simulated power output: For the data training, the hybrid learning algorithm of the ANFIS combines the gradient method with the least squares method to update the parameters [49].In the forward pass of the learning algorithm, consequent parameters are identified by the least squares estimate.In the backward pass, the error signals, which are the derivatives of the squared error with respect to each node output, propagate backward from the output layer to the input layer.In this backward pass, the premise parameters are updated by the gradient descent algorithm.

Multiple Linear Regression Analysis Model
The multiple linear regression analysis model (MLRA) is often used in the prediction being represented by the relationship between inputs and a set of output variables [50].The general equation is as follows: where b is the error vector, which consists of systematic modeling errors and random measurement errors assumed to have a normal distribution and an expected value ( ) 0 E b = .Estimates of the parameter values of 0 , β β are determined by minimizing b.This is simply done using least squares.In the least-squares model, the best-fitting line for the observed outputs is calculated by minimizing the sum of the squares of the vertical deviations from each data point to the line (if a point lies on the fitted line exactly, then its vertical deviation is 0).

Basic Ideas
To explicate the BMA method, let , , , , , denote an ensemble of a prediction obtained from K different models, and Δ be the quantity of interest.In BMA, each ensemble member forecast k f is associated with a condition probability density function (pdf), gk(Δ|fk), which can be interpreted as the condition pdf on Δ at k f , given that k f is the best forecast in the ensemble.The BMA predictive model for dynamic ensemble forecasting can then be expressed as the

Multiple Linear Regression Analysis Model
The multiple linear regression analysis model (MLRA) is often used in the prediction being represented by the relationship between inputs and a set of output variables [50].The general equation is as follows: where b is the error vector, which consists of systematic modeling errors and random measurement errors assumed to have a normal distribution and an expected value E(b) = 0. Estimates of the parameter values of β 0 , β are determined by minimizing b.This is simply done using least squares.
In the least-squares model, the best-fitting line for the observed outputs is calculated by minimizing the sum of the squares of the vertical deviations from each data point to the line (if a point lies on the fitted line exactly, then its vertical deviation is 0).

Basic Ideas
To explicate the BMA method, let f = f 1 , f 2 , . . ., f k , . . ., f K denote an ensemble of a prediction obtained from K different models, and ∆ be the quantity of interest.In BMA, each ensemble member forecast f k is associated with a condition probability density function (pdf), gk(∆|f k), which can be interpreted as the condition pdf on ∆ at f k , given that f k is the best forecast in the ensemble.The BMA predictive model for dynamic ensemble forecasting can then be expressed as the finite mixture model: where w k denotes the posterior probability of forecast k being the best.The w k can be viewed as weights reflecting an individual model's relative contribution to predictive skill over the training period, The value for a k and b k are bias-correction terms that are derived by simple linear regression of ∆ on f k for each of the K ensemble members.The BMA predictive mean can be computed as: This can be viewed as a deterministic forecast whose predictive performance can be compared with the individual forecasts in the ensemble or with the ensemble mean.

The Expectation-Maximization (EM) Algorithm
Successful implementation of the BMA method described in the previous section requires estimates of the weights.There are two widely used approaches for computing the BMA weights, namely the Markov Chain Monte Carlo (MCMC) algorithm and Expectation-Maximization (EM) algorithm.The EM method exhibits many desirable properties as it is relatively easier to implement and is computationally more efficient than the MCMC method [51][52][53].In this study, the EM algorithm was used to identify the BMA parameters.To implement the EM algorithm for the BMA method, the unobserved quantity z kst is adopted, where z kst = 1 if ensemble member k is the best forecast for verification site S and time t, and z kst = 0 otherwise.For each (s, t), only one of {z 1st , . . . ,z Kst } is equal to 1; the others are all zero.
The EM algorithm is iterative and alternates between two steps, the expectation (E) step and maximization (M) step.After initialization of the weights and variances of the individual ensemble members, the EM algorithm alternates iteratively between the E and M step until the convergence is achieved.In the E step, the value of z kst are re-estimated given the current values for the parameters: where g(∆ st f kst , σ j−1 ) is the conditional pdf of ensemble member k with mean f kst and standard deviation σ j−1 at ∆ st , f kst is the kth forecast in the ensemble for location s and time t, and j represents the iteration number.
In the M step, the values of w k and σ 2 are updated using the current estimates of z j kst : Water 2018, 10, 1099 9 of 20 The operated and simulated data were normalized under MATLAB's box-cox function before using the EM algorithm.

Model Performance Metrics
In this study, the Root Relative Mean Square Error (RRSE), Nash-Sutcliffe efficiency coefficient (NSE), and determination coefficient (R 2 ) are used to evaluate the performance of BMA and its three dependent regression models (MLRA, PSO-SVM, and ANFIS).RRSE is calculated as the ratio of the Root Mean Square Error (RMSE) and standard deviation of measured data, is always non-negative, and values equal to 0.0 indicate a perfect fit [54].The value of NSE can oscillate within the interval −∞ ≤ NSE ≤ 1.Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas values <0.0 indicates that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance [55].R 2 ranges from 0 to 1, with higher-values indicating less error variance, and typically values greater than 0.5 are considered acceptable [56].The classification of goodness analysis for NSE and RRME is established in Table 1.The three metrics are defined follows: where N o,t is the observed data set, N s,t is the simulated data set, T is the number of the data set, and N o is the average value of observed value.

Case Study
This study takes Xinanjiang Reservoir as a case study.Xinanjiang Reservoir is located in the upstream of Qiantang River in China, which is the first self-designed and constructed reservoir in China (presented in Figure 3).The basin upstream from the dam site has an area of 10,442 km 2 with a total length of 323 km.Xinanjiang Reservoir is mainly utilized for hydropower generation.The characteristic parameters of the Xinanjiang Reservoir are listed in Table 2, and the conventional operating rules for hydropower generation are shown in Figure 4.The monthly streamflow data collected at the entry of Xinanjiang Reservoir from January 1962 to December 2007 were used as the inflow to the reservoir.To derive the monthly optimal operating rules for Xinanjiang Reservoir, a detailed set of instructions that an operator should follow is described here.
Step 1: Establish an optimal deterministic operation model for Xinanjiang reservoir; Step 2: Set the water level as the decision variable, then use the DP method to solve the model to obtain the monthly power output trajectories (output variables) and reservoir information like the initial water level Z0, inflow Q, maximum water level Za, the incoming hydropower Ef, the storing hydropower Es, and the interaction between incoming and storing hydropower Efs (preselected input variables); Step 3: Select the monthly input variables according to their correlation with the power output using the GRA method from the preselected input variables; Step 4: Derive the individual monthly operating rules using the PSO-LSSVM, ANFIS, and MLRA models.Take different periods for calibration and validation to get the mapping function between input variables and power output as the operating rules with a good regression performance; Step 5: Compute the different weights of the PSO-LSSVM, ANFIS, and MLRA models using the BMA method, and the weighted average results is defined as the final operating rules.To derive the monthly optimal operating rules for Xinanjiang Reservoir, a detailed set of instructions that an operator should follow is described here.
Step 1: Establish an optimal deterministic operation model for Xinanjiang reservoir; Step 2: Set the water level as the decision variable, then use the DP method to solve the model to obtain the monthly power output trajectories (output variables) and reservoir information like the initial water level Z0, inflow Q, maximum water level Za, the incoming hydropower Ef, the storing hydropower Es, and the interaction between incoming and storing hydropower Efs (preselected input variables); Step 3: Select the monthly input variables according to their correlation with the power output using the GRA method from the preselected input variables; Step 4: Derive the individual monthly operating rules using the PSO-LSSVM, ANFIS, and MLRA models.Take different periods for calibration and validation to get the mapping function between input variables and power output as the operating rules with a good regression performance; Step 5: Compute the different weights of the PSO-LSSVM, ANFIS, and MLRA models using the BMA method, and the weighted average results is defined as the final operating rules.To derive the monthly optimal operating rules for Xinanjiang Reservoir, a detailed set of instructions that an operator should follow is described here.
Step 1: Establish an optimal deterministic operation model for Xinanjiang reservoir; Step 2: Set the water level as the decision variable, then use the DP method to solve the model to obtain the monthly power output trajectories (output variables) and reservoir information like the initial water level Z 0 , inflow Q, maximum water level Z a , the incoming hydropower E f , the storing hydropower E s , and the interaction between incoming and storing hydropower E fs (preselected input variables); Step 3: Select the monthly input variables according to their correlation with the power output using the GRA method from the preselected input variables; Step 4: Derive the individual monthly operating rules using the PSO-LSSVM, ANFIS, and MLRA models.Take different periods for calibration and validation to get the mapping function between input variables and power output as the operating rules with a good regression performance; Step 5: Compute the different weights of the PSO-LSSVM, ANFIS, and MLRA models using the BMA method, and the weighted average results is defined as the final operating rules.

Optimal Deterministic Operation for the Hydropower Reservoir
The average water level of 98.4 m was set as the initial operation water level.Figure 5 shows the long-term power output and water level from January 1962 to December 2007, while Tables 3  and 4 represent the monthly results.In this study, one year could be divided into four periods: the flood season (March-July), the transition period from the flood to the non-flood season (August), the non-flood season (September-January), the transition period from the non-flood to the flood season (February).The optimal reservoir operation produced an average power output of 224.71 MW.The power output mainly occurred in the flood season (March to July) with a total power output of 1828.81MW, and a value of 573.22 MW generated during the non-flood season (September to January).The average water level was 106.69 m, and a maximum water level of 107.31 m first appears in August, resulting in enough water stored to generate more power output in the dry period.The minimum water level of 106.14 m in Feb indicated that reservoir could make full use of both natural water resources and reservoir storage to produce power generation in the non-flood season.

Optimal Deterministic Operation for the Hydropower Reservoir
The average water level of 98.4 m was set as the initial operation water level.Figure 5 shows the long-term power output and water level from January 1962 to December 2007, while Tables 3 and 4 represent the monthly results.In this study, one year could be divided into four periods: the flood season (March-July), the transition period from the flood to the non-flood season (August), the nonflood season (September-January), the transition period from the non-flood to the flood season (February).The optimal reservoir operation produced an average power output of 224.71 MW.The power output mainly occurred in the flood season (March to July) with a total power output of 1828.81MW, and a value of 573.22 MW generated during the non-flood season (September to January).The average water level was 106.69 m, and a maximum water level of 107.31 m first appears in August, resulting in enough water stored to generate more power output in the dry period.The minimum water level of 106.14 m in Feb indicated that reservoir could make full use of both natural water resources and reservoir storage to produce power generation in the non-flood season.3. Average monthly power output of hydropower reservoir simulated using the optimal/conventional operation, PSO-LSSVM, ANFIS, MLRA, and BMA models.

Table 3.
Average monthly power output of hydropower reservoir simulated using the optimal/conventional operation, PSO-LSSVM, ANFIS, MLRA, and BMA models.
Average monthly water level of the hydropower reservoir simulated using the optimal/conventional operation, PSO-LSSVM, ANFIS, MLRA, and BMA models.

Monthly Optimal Input Variables Selection Using GRA
The monthly input variables were defined according to the grey relational grade with the mean power output N. The detailed analyses were conducted by season.This can be seen in Table 5: (1) In the flood season (March-July), the monthly average grey relational grade ranged from 0.490 to 0.767 of the input variables of Q, E f , and E fs , which were much higher with a percent of 71.5% to 150.1% than those of the other input variables.Thus, the input variables of Q, E f , and E fs were selected as the optimal factors with the maximum average inflow of 574.32 m 3 /s during this period.The only exception was that all input variables were chosen in May due to both the high-water level and the large inflow.(2) During the transition period between the flood and non-flood season (August/February), the input factors of inflow and reservoir water level had the equivalent influence on the power output.Therefore, all the initial input variables were taken into account here.(3) In the non-flood season (September-January), the hydropower generation mainly relied on the water level and hydropower storage with a decreased average inflow of 107.69 m 3 /s compared to that in the flood season.The input variables of Z 0 , Z a , E f were chosen in this period.

Monthly Operating Rules for Hydropower Reservoir
The monthly average power output and average water level operated based on the power output strategy simulated by the regression models are shown in Tables 3 and 4, respectively.They show that July to October match better than the other months.Though all the modelled monthly water levels are a little bit higher than that of the real optimal ones, no value goes beyond the reservoir upper boundary.
(1) PSO-LSSVM operating rules To derive the SVM operating rules, the PSO was used to obtain optimal values of C and σ, which are listed in Table 6.The monthly simulated power output of PSO-LSSVM generated an average power output of 223.22 MW.The average water level operated by PSO-LSSVM was 107.04 m with an increase of 0.38% compared with that of the optimal results.(2) ANFIS operating rules The ANFIS model produced a second-ranked mean power output of 224.49MW.The average water level simulated by ANFIS was 107.07 m with an increase of 0.39% in comparison to that of the optimal results.
(3) MLRA operating rules The relationship between the monthly power output and the monthly input variables was established for the MLRA model.Table 7 lists the monthly regression coefficients calculated with the "Regression" function in MATLAB.Based on the optimal trajectories, the MLRA generated the maximum average power output of 225.74 MW and highest water level of 107.17 m.Based on the monthly power output trajectories of the three individual reservoir operating rules, the BMA method computed the weights for the PSO-LSSVM, ANFIS, and MLRA models.The weights reflect the performance of ensemble models.As can be seen from Table 8, MLRA was the worst performing among the ensemble models for simulating the operating rules, while the other two models trajectory provided the targeted water level operating for the PSO-LSSVM, ANFIS, MLRA, and BMA models.The RRSE and NSE of all simulating models calculated for the long-term water level are also shown in Table 2.All the regression models were within the very good range (NSE > 0.6, RRSE < 0.5).The BMA performed a first rank with the maximum NSE of 0.724 and the minimum RRSE of 0.185 during the calibration period and the maximum NSE of 0.956 and the minimum RRSE of 0.179 during the validation period.Therefore, the BMA model could achieve a relative optimum.

Comparison with the Conventional Operating Rules
Similar to the optimal operation procedures, the monthly input variables were determined according to the grey relational grade with the mean power output under the conventional operation.Seen from Table 11, the simulating power output was mainly dominated by the water level during the whole year.There was an exception in July due to its large inflow.The conventional operating rules produced a reduced mean power output of 192.63 MW compared to the optimal and simulated ones, the same situation to the average water level of 100.55 m.This states that the optimal operating rules can make better use of the natural inflow and improve the power generation efficiency compared with the conventional rules.

Comparison with the Conventional Operating Rules
Similar to the optimal operation procedures, the monthly input variables were determined according to the grey relational grade with the mean power output under the conventional operation.Seen from Table 11, the simulating power output was mainly dominated by the water level during the whole year.There was an exception in July due to its large inflow.The conventional operating rules produced a reduced mean power output of 192.63 MW compared to the optimal and simulated ones, the same situation to the average water level of 100.55 m.This states that the optimal operating rules can make better use of the natural inflow and improve the power generation efficiency compared with the conventional rules.

Conclusions
Various regression models have been applied to the derivation of operating rules.It is necessary to analyze and evaluate the model selecting uncertainty involved in reservoir operating rules.Moreover, selecting the optimal input variables from a large number of candidates to characterize an output variable can lead to a more accurate operation simulation.For efficient hydropower generation, this study coupled the GRA and BMA methods to derive monthly optimal operating rules for reservoir power generation.The primary processes were as follows: (1) an optimal deterministic operation model of reservoir power generation was formulated and solved; (2) based on the monthly optimal deterministic operation strategies, the monthly input variables were selected according to their correlation with the power output using the GRA method; (3) PSO-LSSVM, ANFIS, and MLRA models were used to derive the individual monthly operating rules; and (4) the BMA was applied to determine the final reservoir operating rules by analyzing the uncertainty of selecting individual models with different weights.The Xinanjiang Reservoir in China was taken as a case in this study, which shows: (1) Inflow, storage, and their formative input variables inconsistently acted in simulating operating rules as their monthly grey relational grade variously ranged during different periods.
(2) The MLRA model performed worst in simulating the operating rules, while the other two models showed similar performance because MLRA model generated the lowest weights compared to the other two models for every month.(3) The BMA method outperformed among all operating rules as it considered model selecting uncertainty and had the ability to improve reservoir operations: (a) the average power output and water level under BMA were approximate to that of the optimal one, and (b) the R 2 of the BMA was greater than that given by any of the individual operating rules, where R 2 reflected the goodness of fit to the optimal trajectory.(4) The combination of the two methods could achieve larger hydropower generation and make better use of natural inflows in comparison to the conventional operation as it produced a larger mean power output than the conventional operating rules.
However, this study only considered one single reservoir; further work will focus on extending the operating rules of a single reservoir to cascade reservoirs or mixed multi-reservoir systems.

Figure 1 .
Figure 1.Flowchart for the derivation of GRA-BMA operating rules.

Figure 1 .
Figure 1.Flowchart for the derivation of GRA-BMA operating rules.

Figure 4 .
Figure 4. Conventional operating rules for hydropower generation of Xinanjiang Reservoir.

Figure 4 .
Figure 4. Conventional operating rules for hydropower generation of Xinanjiang Reservoir.

Figure 4 .
Figure 4. Conventional operating rules for hydropower generation of Xinanjiang Reservoir.

Figure 5 .
Figure 5.The hydropower reservoir operation results of the optimal/conventional operation from January 1962 to December 2007.

Figure 5 .
Figure 5.The hydropower reservoir operation results of the optimal/conventional operation from January 1962 to December 2007.

Figure 6 .
Figure 6.Scatter plots of the monthly power output given by PSO-LSSVM, ANFIS, MLRA, and BMA under the optimal trajectory.

Figure 6 .
Figure 6.Scatter plots of the monthly power output given by PSO-LSSVM, ANFIS, MLRA, and BMA under the optimal trajectory.
ensemble members can be approximated by a normal distribution centered at a linear function of the origin forecast, a k + b k f k , with standard deviation σ k

Table 1 .
Classification of goodness of fit for NSE and RRSE.

Table 2 .
The characteristic parameters of the Xinanjiang hydropower station.

Table 5 .
The monthly input variables defined based on the optimal trajectory.

Table 6 .
The monthly optimal values of C and σ of LSSVM using PSO algorithm.

Table 7 .
The monthly regression coefficients for the MLRA model.

Table 9 .
The results of NSE and RRSE for the monthly power output simulated by PSO-LSSVM, ANFIS, MLRA, and BMA models.

Table 10 .
The RRSE and NSE of the long-term power output/water level simulated by PSO-LSSVM, ANFIS, MLRA, and BMA models.

Table 11 .
The monthly input variables defined based on the conventional trajectory.

Table 11 .
The monthly input variables defined based on the conventional trajectory.