Hybridizing ANN-NSGA-II model with genetic programming method for reservoir operation rule curve determination (Case study Zayandehroud dam reservoir)

One of the most important and effective works of water resource planning and management is determining the specific, applicable, regulated operating policies of the Zayandehroud dam reservoir, as a case study, in which it should be user-friendly and straightforward for the operator. For this purpose, different methods have been proposed in which each of them has its limitations. Due to the unique capabilities of the genetic programming (GP) model, here, this method is used to determine the operating rule curves and policies of the dam reservoir. For this purpose, here, two cases are proposed in which, in the first case, each month is individually simulated and modeled. However, in the second case, all months are simulated simultaneously. A second case is proposed here to determine simple and more applicable operation rule curves. In addition, two approaches are suggested for each case in which in the first approach, the influential input variables are selected by presenting the hybrid method. In the proposed hybrid method, the artificial neural network (ANN) model is equipped with non-dominated sorting genetic (NSGA-II) algorithm leading to a hybrid method named the ANN-NSGA-II method. However, in the second approach, the influential input variables are selected automatically using the GP method. Here, the hybrid method is proposed and used to overcome the limitations of existing usual method. In other words, it is proposed to reduce the number of influential input variables of data-driven methods and select the effective ones. The obtained results of all proposed cases and approaches are presented and compared with the standard operation policy method, stochastic dynamic programming, ANN model, and NLP method. Comparison of the results shows the acceptable performance of the proposed cases and approaches. In other words, the best-obtained values of (stability index) SI index and water deficit (objective function value) are 49.3% and 32, respectively.


Introduction
Nowadays, the water resources limitation becomes a sever problem in many areas which is restricted the area's developments. Therefore, water resource management is an essential challenge for these areas. Generally, water storage of dam reservoirs is one the most critical water resource. However, determining the optimal operation policies is a complex decision process in which the optimal amount of water release from a reservoir is determined based on the reservoir storage capacity, water demand, and the inflow into the reservoir. Generally, to optimally operate the water stored in the dams' reservoirs, it is necessary to determine the reservoir operation policies presenting as rule curves or tables leading to determine the optimal water release or the storage volume values of a reservoir at different operational time period.
Due to this fact, these days, many scientific and applicable types of research have been done in the research field of operation rule curve and policies determinations. For example, Karamouz and Houk (1982) used the dynamic programming (DP) method, regression analysis, and a simulation model to obtain operating rules and policies of the dam reservoir considering 48 states. In this study, the annual operating rules were determined considering 12 states, and the monthly operating rules were determined considering 36 states. The results showed that this proposed method was a very effective and simple method for reservoir rule curve determination. Karamouz and Houk (1987) compared the performance of two methods, which means the DP and stochastic dynamic programming (SDP), to determine the reservoir operation rule curves. Comparison of the results showed that the performance of the DP method was better for determining the operation rules of large-scale reservoir and the performance of the SDP method was better for small cases. Rossi et al. (1999) used the discrete differential dynamic programming (DDDP) method with the objective function of minimizing the deficiency of water releases and demand to optimize the water release values from the reservoir. In addition, multiple linear regression method and artificial neural network (ANN) model were also used to determine the reservoir operation policies. Here, in these methods, water storage values at the start of the operation time period, water inflow and water releases values were considered as state variables. The results showed that using the ANN method was improved the efficiency of the method in comparison with the regression method. Mousavi et al. (2004) proposed a new method named fuzzy stochastic dynamic programming (FSDP) by combining the SDP and fuzzy logic methods. In this proposed method, the uncertainty of the hydrological variables and the inaccuracy of the results of the discrete variables were investigated. In this method, the transition probability matrix was obtained using the fuzzy logic and the fuzzy Markov chain. The results of the proposed hybrid FSDP method were compared here with the results of the SDP method, in which the results of the FSDP method were better than the SDP method. Chang et al. (2005) used the genetic algorithm (GA) and fuzzy logic method to solve this problem. In other words, they used the adaptive network-based fuzzy inference system (ANFIS) model to determine the operation rule curves of the reservoir. The results showed that using artificial intelligence improved the system performance and the ANFIS model was a capable model for the learning process of different data and information resources. Kerachian and Karamuoz (2006) used water quality simulation models and the deterministic and stochastic models to determine the reservoir operation rule curves and policies. In this research, Nash's bargaining theory was used to resolve the stakeholder's conflicts. The simulation model was used here to obtain the thermal stratification of the reservoir and the quality of the water release from the reservoir. In addition, the varying length genetic algorithm (VLGA) was used in this study, in which the computational complexity is less than the SDP model. The results of this study showed that the proposed model was a capable model to reduce the salinity of the water release from the reservoir and water storage in the reservoir. Paredes and Lund (2006) used a linear programming (LP) method to determine the reservoir operation rule curves considering the water quality parameters during drought and wet seasons. Chang and Chang (2009) used the non-dominated sorting genetic algorithm (NSGA-II) algorithm to determine the operation rule curves of the multi-objective multi-reservoir system. The results of this research indicated the capabilities of this algorithm for determining the operation rule curves of the multi-objective multi-reservoir system. Karamouz et al. (2009) proposed a new quantitative-qualitative model for the stochastic operation of the dam reservoir considering the uncertainty of water inflow into the reservoir. In this study, the operation rule curves were determined and the water allocation problem was solved to increase the reliability of water quality at the downstream and upstream of the reservoir. For this purpose, a method, named the Bayesian stochastic genetic algorithm (BSGA) model, was used to determine the reservoir operation rule curve considering the uncertainties of water inflow values. The obtained results showed the capability of this proposed method. Malekmohammadi et al. (2009) used a model named, varying chromosome length genetic algorithm (VLGA-II) and Bayesian method, to determine the reservoir operation rule curves. The proposed approach was used to minimize the long-term water scarcity damages and short-term flood damage of multi-reservoir systems. Here, the obtained operation rule curves using the Bayesian network were compared with the fuzzy linear regression model. Comparison of the results showed that the loss values were reduced by using the Bayesian theory. In addition, the obtained results of this method were more accurate than other methods. Jothiprakash et al. (2011) used GA and SDP methods to determine the operation rule curves of a multi-reservoir system. Here, minimizing the deviation of water demand was considered as objective function. Comparison of the results showed that the performance of the GA was better than the SDP method. In addition, the water reseals values obtained with the SDP model were not able to supply the water demand for most months. Castelletti et al. (2013) used a reinforcement learning algorithm, and a method called fitted Q iteration (FQI) to determine the operation rule curves of a multiobjective reservoir. The results showed that using this approach increased the water quality downstream and in the dam reservoir. In addition, the water allocated for agriculture purposes was extremely reduced. Fu et al. (2013) used a method called factorial-based fuzzy stochastic dynamic programming to solve the reservoir operation optimization problem. Here, the proposed algorithm was incorporated stochastic, fuzzy processes and their relationship. This algorithm was a combination of the SDP method, fuzzy Markov chain, and factorial analysis techniques. This model was an applicable model to determine the optimal water release policies from reservoir at uncertain conditions. Boluri-Yazdeli et al. (2014) used the standard operating policy (SOP), SDP, linear decision rule (LDR), and nonlinear decision rule (NLDR) methods to determine the operation rule curves of the dam reservoir. In this study, a useful multi-criteria decision-making method, which means the ELECTRE-I method, was used to determine the operation rule curve. To evaluate the used methods, the reliability and residency indexes and the system efficiency were calculated. This study showed that the residency of LDR and that of NLDR methods were the lowest values and reliability values were highest and, overall, using the NLDR method led to better results. Taghian et al. (2014) used GA and LP methods to determine the operation rule curves of multi-objective multireservoir system for water resources restriction conditions. Here, the objective function of the LP model was to minimize the cost function. The results showed that the search space size of the GA was reduced by using the LP model. In addition, this model was a capable model to model the droughts periods over several years. Akbari Alashti et al. (2014) used the nonlinear programming (NLP) method, GA, and fixed-length gene genetic programming (FLGGP) for the online operation of the three-reservoir system. Comparison of the results showed the better performance of the FLGGP method for determining the reservoir operation rule curves. Li et al. (2014) used the genetic programming (GP) method to determine the operation rule curve of a multi-reservoir system. For this purpose, at first, the optimal operation policies were determined using the DP method. Comparison of the results showed that the GP results were better than the ANN results, and energy production and system reliability were increased by using the GP method. Ashofteh et al. (2015) used GP to determine the reservoir operation rule curve for normal and climate change conditions. Here, water inflow, reservoir storage volume, and downstream water demand values were considered as GP model input variables. In addition, maximization of system reliability and minimization of system vulnerability were considered as objective functions of the problem. The results showed that the GP results were improved 29% up to 32% by considering climate change conditions. Zhang et al. (2015) used the Bayesian model averaging (BMA) to determine the operation rule curves of the dam reservoir considering uncertainties of the problem. The proposed model consisted of three steps as follows: (1) deterministic optimization using the DDP2 and (2) using the least square support vector machine (LS-SVM), twodimensional surface model (SURF) and piecewise linear regression (PL-REG) models to determine the operation rule curve, and (3) combining the three obtained operation rule curves. The results showed that the water release values were increased by using the DDP2 method in comparison with the results of the PL-REG, LS-SVM, and the BMA models. Spiliotis et al. (2016) used the particle swarm optimization (PSO) algorithm to determine the hedging rule of a multi-reservoir system. Two models were proposed here. The first model consisted of a hedging rule curve that only included drought and no drought conditions, and a specific bound between the two conditions was also considered. However, in the second model, two hedging rule curves were defined, including no drought, alert, and emergency conditions. The starting point for solving the problem was risk analysis to determine the boundary that it limited the water release from the reservoir. The proposed mechanism was used to reduce the problem dimension. Sharma et al. (2016) used the SDP method to determine the optimal operation policies for dam reservoirs. In this research, the water storage volumes of the reservoir at the start of the operation time and the water inflow into the reservoir were considered as the state variables. In addition, minimizing the difference between the water storage volumes and the water release values and their expected values were considered as the objective function of the problem. In addition, the effect of state variable discretization for determining the optimal reservoir operation policies was also investigated. The results showed that if the unequal discretization intervals were considered for water storage volume variables, the objective function values were improved 8% to 58% compared to similar discretization intervals method. Finally, the water storage volumes at the end of each operation time were determined for different conditions. Haguma et al. (2018) used four methods called counting, ordinary leastsquares regression, robust linear model regression, and multivariate conditional distribution to determine the transition probability matrix and finally to study the effect on the short-term operation of the multi-reservoir system. In this research, two approaches of the equal interval and the equal number of discretization were used to discretize the variables. Comparison of the results showed that the discretization approaches had the highest effect on the results of the regression and counting methods, and the least effect on the multivariate conditional distribution method. Saadat and Asghari (2017) proposed a new method named reliability improved stochastic dynamic programming (RISDP) method to determine the reservoir operation rule curves. In proposed method, the infeasible solution was not created compared with the usual SDP method. Here, the concept of reliability was used to increase the water release values from the reservoir to satisfy water demand. The results showed that using this algorithm was improved the results and avoided creating an infeasible solutions. Lei et al. (2018) proposed a new and more accurate method for determining the transition probability matrix for reservoir operation policies determination using the SDP method. In this study, the Copula function was used to determine the transition probability matrix, and the function value was determined at the last operation period using an iterative method. This method was studied for a dam reservoir in China, and its results were compared with the conventional and usual SDP method. Saadat and Asghari (2018) determined the optimal water storage volumes of the reservoir using SDP and a new method named cooperating SDP and nonlinear programming (CSDP). The acceptable results were obtained for this model by using a smaller number of state variables. The results showed that the reliability of the CSDP model was better than that of the SDP model. In addition, the distance-based interpolation formula (DBIF) was used to evaluate the efficiency of different methods proposed for determining the reservoir operation policies. The results showed that the performance of the CSDP model was better than the multiple linear regression (MLR) model. Ashofteh et al. (2019) used bi-objective genetic programming (BO-GP) algorithm to optimize the reservoir operation rules and policies. Here, the BO-GP algorithm was able to determine the optimal reservoir operation policies by calculating the Pareto front. The results showed the acceptable performance of the BO-GP algorithm to solve the multi-objective reservoir operation optimization problem. In addition, the system vulnerability and reliability were 16% up to 41% and 46% up to 78%, respectively.
Investigating the mentioned researches indicates that proposing new methods or modifying the existing techniques to determine the reservoir operation rule curves and policies is one of the most significant and attractive fields for water resources researchers and operators. Due to this fact, the applicable and straightforward form of operation rule curves is requested by the operator, planner, and manager for sustainable water resource management. Due to the unique characteristics of the GP method, therefore, the operation rule curves of the Zayandehroud dam reservoir are determined here using the GP method by proposing two cases. In the first case, the operation rule curve of each month is obtained individually. However, in the second case, the operation rule curves of all months are obtained simultaneously. These proposed cases are one of the innovations of this research, leading to more applicable operation rule curves. In addition, here, two different approaches are also presented for each case. In the first approach, the influential input variables are selected by proposing the hybrid ANN-NSGA-II method, combining the ANN model with NSGA-II (non-dominated sorting genetic) algorithm. However, in the second approach, the influential input variables are selected automatically using the GP method. Reviewing the researches shows that some hydrological variables such as water discharge and various conditions of the reservoir such as water storage volume and water demand at the downstream are influential and essential variables and factors which affect the reservoir operation policies for water supply. Therefore, here, the water inflow, precipitation, demand, and storage volume values with different time lags (1-12 months) are initially considered here as input variables. In addition, another input variable means time index T related to each operation time is also considered to improve the accuracy of obtained results. It should be noted that one of the most crucial challenges of data-driven methods such as GP is selecting the useful and practical input data set. Therefore, different techniques such as experimental and trial and error methods have been proposed to solve this problem. Each of these methods has its limitations and disadvantages. Therefore, here, the hybrid ANN-NSGA-II method is proposed and used to reduce the number of influential input variables and select the useful ones. This fact is another innovation of this research, which is fully presented in Sect. ''Proposed method.'' As a case study, the operation policies of the Zayandehroud dam reservoir are predicted using proposed methods and compared with the obtained results of other methods. For comparison propose, all obtained results are compared here with the SOP (standard operation policy) method, SDP, ANN model, and NLP methods using different indexes such as SI (sustainability index). Comparison of the results shows the acceptable performance of the proposed cases and approaches highlighted in Sect. ''Results and discussion.'' Therefore, the structure of the present research is as follows. By reviewing the literature and highlighted the research innovations in the introduction section, in the second section, a brief description of the used methods is presented. The data and information of the study area (Zayandehroud dam reservoir) are shown in the third section. In the fourth section, the proposed methods are presented. The results of the proposed methods are presented and compared with other methods such as ANN in the fifth section. Finally, in the sixth section, some concluding remarks are given.

Basic concept of the used method
In this section, a brief description of the used method is presented. Here, the GP method is used to determine the operation rule curves of the Zayandehroud dam reservoir, and therefore, it is briefly described. By combining the ANN model and NSGA-II algorithm, a method is proposed for selecting influential input variable called the ANN-NSGA-II method. Therefore, in this section, the ANN model and NSGA-II algorithm are briefly described. In addition, the operation rule curves of the Zayandehroud dam reservoir are determined using SOP, SDP, ANN, and NLP methods for comparison purposes, and therefore, these methods are also briefly described.
It is worth noting that these methods are usual traditional methods which are fully described in many types of research. Therefore, details of these methods are not fully presented in this paper. In this section, brief descriptions of these methods are only given.

Genetic programming (GP)
In general, based on the concept of the genetic algorithm (GA), the GP method was proposed, which is inspired using the naturally evolutionary rule. Generally, finding a direct relationship between physical parameters is a unique future of the GP method in comparison with other datadriven method leading to an applicable approach (Shiri and Kishi 2011). Furthermore, influential input variables can be easily selected using this method (Naseri et al. 2011). The basic steps of using the GP method to solve problems are presented as follows: (1) The initial population is selected randomly.
(2) Best chromosomes are selected and repeated for next generation. The roulette-wheel, tournament, and ranking methods are more common selecting methods.
(3) To start the evolution process of generations, better chromosomes are used, considering crossover and mutations operators with predefined probabilities. It should be noted that a function and terminal can be replaced with another function and terminal, respectively, in this process. Generally, by using the crossover operator, some parts of the chromosome are replaced with the same parts of the other chromosomes. Furthermore, using the mutation operator, a function or number can be replaced with another function or number, respectively. (4) Among finding solution processes, better chromosomes are replaced with worst chromosomes, and best chromosomes remain unchanged. (5) These steps are continued until the best solution is obtained (Johari et al. 2006).
It is worth noting that the terminal consisting of random numbers and variables, the mathematical functions, and appropriate fitness function should be set before the GP process is started.
The general idea of the GP is to use parse trees as chromosomes which can be capture expressions in a given formal syntax. Depending on the problem, and the users' perceptions of what the solutions should look like, this can be the syntax of arithmetic expressions, formulas in firstorder predicate logic, or code written in a programming language based on the parameters set.
Due to the unique characteristics of the GP, this model was highly used for a variety of machine learning tasks in different fields. A review of the GP model and its applications in the field of civil engineering was presented by Zhang et al. (2021) over the last decade. In addition, some modified forms of GA were also proposed in different fields of engineering problems to overcome the limitations, such as researches of Wu et al. (2011), Pawlak et al. (2015, and Bourisli et al. (2018).

Non-dominated sorting genetic (NSGA-II) algorithm
Based on the Darwinian evolution theory and the natural evolution process, a new algorithm named GA was proposed. In this theory, only the survivals of creatures that are more consistent with their surroundings and more graceful remain. In this algorithm, each decision variables values represent by a gene, and the combining of genes leads to a chromosome. Generally, the GA process to solve an optimization problem can be summarized as follows: 1-The initial population of chromosomes is selected. 2-The fitness function values of each chromosome are determined. 3-Chromosomes with better fitness function values are selected for the next generations. 4-The crossover and mutation operators for subsequent population are applied. 5-The fitness function values of the new population are determined. 6-Steps 2 to 5 are continued until the stopping criteria are reached.
Generally, using the concept of the GA, the NSGA-II algorithm was proposed. In this algorithm, based on the number of dominated solutions in comparison with other existing solutions, a rank has been assigned to each solution. In other words, a rank one has been assigned to the solutions that have not been dominated, representing the Pareto front rank one. In continue, the Pareto front rank one solutions are removed from search space, and remained solutions are compared with each other leading to Pareto front rank two. In other words, rank of two has been assigned to the solutions that have not been dominated. Finally, this process is continued until all solutions are classified into different Pareto front with other ranks.
The crowding distance concept is generally used by the NSGA-II algorithm to decide which points of a particular Pareto front will be removed or remained, presenting by Eqs. (1), (2), and (3). Therefore, a D value is calculated, and the points with smaller D values are eliminated for all points of the particular front. The more considerable D value represents that, removing these points might decrease the greater variety of the search space. Generally, the NSGA-II algorithm process to solve a multi-objective optimization problem can be summarized as follows: 1-The initial generation,P t , is generated in which it contains N members. 2-The objective function values of each chromosome are calculated. 3-The parent generation is sorted, and each solution is ranked based on the non-dominant sorting approach. 4-By using crossover and mutation operators for the initial generation, the new generation, Q t , is produced. 5-By combining P t and Q t generations, the new generation, R t , is produced in which it contains 2 N number. 6-Sorting the new generation,R t , is done based on the non-stack classification method, and therefore, Pareto front with a different rank is produced. 7-Based on the crowding distance concept, N members are selected to create P tþ1 generation from R t generation in which it has 2 N members. 8-The created generation, P tþ1 , is sorted, and the crossover and mutation operators are used to produce the new generation, Q tþ1 . 9-This process is repeated until the best solution is obtained (Sepahvand et al. 2019).

Artificial neural network (ANN) model
Based on the human brain neural network, an intelligent computational system was proposed, which called the ANN model. This model comprises a combination of information processing units called neurons in which one or more inputs are received, and the outputs are determined using the predetermined non-linear function (Raman and Sunilkumar 1995). Nowadays, the ANN model solves many problems due to unique features (Ferench et al. 1992).
To create the structure of an ANN model, the number of neurons, input, hidden and output layers, and weights should be defined. Different methods are proposed and used to determine the proper values of these parameters. Generally, these parameters can be determined using experimental and try and error methods to avoid under and, or over fitting problems. In addition, many different ANN models have been proposed by combining these parameters. One of the most famous and straightforward ANN models is feed-forward neural network (FFNN) model, which is used to predict time series variables in different fields such as water engineering using non-linear functions for hidden layers. Nowadays, various forms of ANN models have also been proposed for this purpose (Nourani et al. 2019;Sharghi et al. 2018). Furthermore, various methods have been generally proposed to reduce the error of the ANN model. Backpropagation network (BPN) is one of the most usual methods in which it commonly used the steepest gradient method to improve the weight of connected neurons. Furthermore, the BPN network is a capable model for enhancing the interaction of processor elements by adding hidden layers (Mishra and Desai 2006).
Another essential parameter of the ANN model is the transition function in which different function such as linear, nonlinear, sigmoid, sigmoid tangent, and so on have been proposed and used for each layer of the ANN model. The proper function can be defined using experimental and trial and error methods. Based on the previous researches, here, the most usual and simple form of the ANN means feed-forward backpropagation network is used for prediction process using the sigmoid function 1 1þe Àx . It is worth noting that all data should be normalized to standardize different values of input and output variables using Eq. (4).
where X = original data; X n = normalized data; X min = minimum value of data; and X max = maximum value of data.

Stochastic dynamic programming (SDP) method
SDP method is the most commonly usual method in which it is completely presented in many types of research, and therefore, a brief description of this method is presented in this section. SDP is a modified form of the dynamic programming (DP) method which is used to solve the stochastic problem. In this method, the state ((s n )) and decision (d n ) variables should be defined and discretized. This method is used to solve serial nature problems such as reservoir operation optimization problem. For this purpose, the problem can be divided into N sub-problem leading to the N stage, and each stage, n (n = 1,…,N), can be solved separately due to Bellman principles. For the reservoir operation optimization problem, each operation period (t) can be considered as a stage. In addition, the state variables are the water storage values at the start of operation period (s(t)). Furthermore, the water storage values at the end of operation period (s(t ? 1)) or water releases for the reservoir at operation period t (R(t)) can be considered as decision variables. Due to the stochastic nature of water inflow into the dam reservoir, in the SDP method, the transition probability matrix should be computed by discretizing the water inflow bounds. By defining the recurrence relationship (function), r, the best decision values for each state variable of each stage can be computed, and this process is continued for the next state until all N stages are covered. This relationship can be defined based on the objective function of the problem, previous decision, and transition probability matrix values. It is worth noting that the state variable of next stage (n ? 1) can be determined using the transition function (s t?1 = t(s t , d t )). Figure 1 shows the SDP process for solving the optimization problem schematically.
It should be noted that the primary and usual form of the SDP method has some limitations and problems. Therefore, here, a method is proposed and used for the SDP method to solve a fundamental problem of this method. In other words, the feasible operation policies cannot be determined for some cases using the SDP method. This problem is solved by considering the maximum or minimum allowable water storage volume values and recomputing the water release values when the dam reservoir operation problem is simulated.

Standard operation policy (SOP) method
The SOP method is one of the usual traditional methods for determining the operation policies of the dam reservoir. In this method, the water releases values from the reservoir proportionally related to the water inflow availability. In other words, it is assumed both past and future conditions, and therefore the water inflows will be available at any time. In this method, if the dam reservoir is empty, the water release value is equal to zero. In addition, if the water storage values at the operation period are greater than the water demand, the water reseal values from the reservoir are similar to the water demand values and the overflow values stored in the reservoir. Furthermore, when the reservoir water storage is not fully capable of storing, some percentage of the water demand is not satisfied (Emadi et al. 2016). Therefore, the water release values and the water storage volumes are calculated based on the following equations: where RðtÞ = water release from the reservoir at operation period t; SðtÞ = water storage volume at the start of operation period t; DðtÞ = water demand at the operation period t; IðtÞ = water inflow into the reservoir at the operation period t; and S max = maximum value of water storage volume.

Fig. 1 SDP process to solve problem
Hybridizing ANN-NSGA-II model with genetic programming method for reservoir operation rule… 3 Case study To investigate and evaluate the performance of the different proposed methods for reservoir rule curve determinations, here, the Zayandehroud dam reservoir is considered as a case study. Generally, the Zayandehroud catchment's area is located between latitudes from 31 15 0 to 33 45 0 north and longitudes from 50 02 0 to 53 45 0 east. The annual average value of rainfall is 140 mm (mm) in which it is varied from at least 50 mm (mm) in the east up to a maximum of 1500 mm in the west (Babaei et al. 2019). This catchments area is 1.6% of the total area of Iran, with a value of 26,972 Km 2 . Figure 2 shows the location of several hydrometric stations of the Zayandehroud catchment (Safavi et al. 2015). The north and south limitations of this area are the Salt Lake catchment and Abarqoo Desert basin, respectively. In addition, the east and west limitations of this area are Dagh-Sorkh Playa and Siyah-Kooh mountains and the Karun River basin, respectively (Babaei et al. 2019).
Zayandehroud dam catchment area is about 4120 Km 2 in which the area of the main branch of the Zayandehroud catchment located at the Ghaleh-Shahrokh station is 1500 Km 2 . In addition, the area of the Pladsjan River catchment located at the Skandari station is 1,600 Km 2 , and the area of the Samandang River catchment located at the Mendarjan station is 227 km 2 . Finally, the area of catchment located after underground hydrometric stations or directly drained into the dam's lake is 793 km 2 (Safavi et al. 2015).
In this research, the observation data from 1976 to 2013 have been used containing the data of Ghaleh-Shahrokh, Skandari, and Mendarjan stations and the natural water inflow (without considering the flow of transfer tunnels) and precipitation. Here, the average annual total precipitation is about 140 mm in which it is varied from 1500 mm in the west to 50 mm in the east (Saadat and Asghari 2017).
It should be noted that the maximum and minimum water storage volumes of the reservoir are 1447.4 cubic per meter (MCM) and 115.2 MCM, respectively. In addition, the maximum and minimum water releases values from the reservoir are 500 and zero MCM, respectively. The observed water inflow values into the Zayandehroud dam reservoir are presented in Fig. 3.

Proposed method
In this research, different methods are proposed to determine the reservoir operation policies of the Zayandehroud dam reservoir for water supply. For this purpose, at first, the influential input variables and essential parameters should be defined for modeling. Reviewing the researches shows that some hydrological variables such as water discharge and various conditions of the reservoir such as water storage volume and water demand at the downstream are influential and essential variables and factors which are affected the reservoir operation policies for water supply. Therefore, for data-driven models, the water inflow, precipitation, demand, and storage volume values with different time lags (1 to 12 months) are initially considered here as input data set. In addition, another input data set, Fig. 2 Location of the Zayandehroud Basin and its sub-basins (Safavi et al. 2015) which means time index T, is also considered to improve the accuracy of the obtained results. This parameter is related to each time of reservoir operation policy. In addition, the water storage volume at the end of the operation period of the target month is considered as the output data set.
In this research, at first, different optimization methods such as SDP and NLP methods are used to determine the optimal operation rule curves of the Zayandehroud dam reservoir. Here, the best-obtained results are selected by evaluating these methods, and finally, these results are used to determine reservoir operation rule curves using GP and ANN methods. It should be noted that the formula is the main output results of the GP method which can be used as reservoir operation policies. Therefore, in the GP method, a relationship is generally superior to select as the best formula in which it involves all or more parameters. In addition, to select the best formula, the RMSE and R 2 values should also be used. Therefore, for the proposed GP model, the formula with the lowest error value and including the maximum number of independent parameters has been selected. In this paper, the equations of section five are the final output of the GP method using proposed cases in which they can be used as reservoir operation rule curves.
In this research, two ceases are proposed for each GP and ANN method. In the first case, each month is modeled separately, and in the second case, all months are modeled simultaneously. Selecting the influential input data set for data-driven methods is one of the most essential challenges of these methods in which different methods have been proposed to solve it. Experimental and try and error methods are the most usual methods used for this purpose. However, here, two approaches are proposed for each case to select the influential input data set. In the first proposed case, the influential input variables are selected by proposing the hybrid ANN-NSGA-II method. However, in the second one, the influential input variables of the GP methods are automatically selected the using specific future of this method. Figure 4 schematically shows the different defined processes of this research.
Generally, considering each of the input variables separately and using the try and error method to investigate its effect on the output results is a usual and formal method to select input data sets for data-driven models. However, a hybrid ANN-NSGA-II method is proposed here for this purpose which is fully described as follows. In the proposed method, the independent physical variables X 1 to X n are assumed and used to determine the dependent variables as output results to select the characteristics and features of variables. The goal of this process is to define a subset of the feature. This feature should be defined as small as possible to predict the new output O correctly. Generally, usual definite approach cannot be proposed for this problem. However, different optimization algorithms can be used for it. This problem is a multi-objective optimization problem in which it can be defined as follows: wherex i j j= the number of chosen features,x i = the chosen inputs, O = the output vector, andf = the objective function of the model. It is worth noting the objective function is generally obtained based on the chosen inputs.
In this research, the NSGA-II algorithm is used to solve this optimization problem. This algorithm is combined here Hybridizing ANN-NSGA-II model with genetic programming method for reservoir operation rule… with the ANN model, and therefore, the hybrid ANN-NSGA-II method is proposed. In this proposed method, a relationship between input and output variables can be found using the ANN model.
In the research, the operation policies of the Zayandehroud dam reservoir are determined using the SDP method. For this purpose, all state variables (water storage volumes at the beginning of the operation period) are discretized into different interval values. The computation process requires a specific number defined as an interval index which is usually the average value of the interval bounds. A review of this research field shows that several approaches have been proposed for discretizing the state variables. These approaches are classified as (1) classic method, (2) Savarenskiy method, (3) Moran method, (4) equal interval values method, and (5) an equal number of data method (Karamouz et al. 2003). Using the equal interval values method leads to some intervals with no data, and therefore, the computational process is disrupted. Due to this fact, in this research, the equal number of data method is used to discrete the state variables. By discretizing the state variables, the Markov matrices determination is another critical step. In this research, an efficient method is used to calculate the probabilities of Markov matrices values. For this purpose, at first, the monthly water inflow value is ascended. Then different intervals (such as two, three, four, or five intervals) are determined so that the number of data in each interval is the same (Saadat and Asghari 2017). Furthermore, here, the operation policies of the Zayandehroud dam reservoir are determined using the NLP method. For this purpose, two linear and nonlinear optimization models are proposed for determining the reservoir operation policies using the historical and observational data (current conditions). Decision-making problems, such as reservoir operation optimization problem, can be defined as a mathematical model. This model can be solved using different methods. For this purpose, it is necessary to determine the decision variables, constraints, and the objective function of the problem.
Generally, the objective function of the reservoir operation problem can be defined as minimizing the cost or loss or maximizing the benefit. In this research, minimizing the summation of water release losses in comparison with the downstream water demand is defined as the objective function of the problem (Soghrati and Moeini 2020): where Z = the objective function of the problem, n = the total number of the operation time period, RðtÞ = water release from the reservoir at operation time period t, DðtÞ = water demand at the operation time period t, and D max = maximum water demand.
Here, the water release from the reservoir at each operation time period is considered as the decision variable of the problem. In addition, all constraints should be defined to define the mathematical optimization model completely. Generally, the continuity equation is the most important constraint of this problem which is defined as: S(t þ 1) = S(t) þ IðtÞ À RðtÞ À LðtÞ ð 8Þ where Sðt þ 1Þ = reservoir storage volume at the end of operation time t (start of operation time t ? 1), SðtÞ = reservoir storage volume at the beginning of operation time t, RðtÞ = water release from the reservoir at operation time t, IðtÞ = water inflow at operation time t, and LðtÞ = water loses at operation time t including infiltration from dam reservoir and evaporation. Other constraints of this problem are water storage and water release bounds which are defined as: where S max = maximum water storage volume of the reservoir, S min = minimum water storage volume of the reservoir, R max = maximum water releases from the reservoir, R min = minimum water releases from the reservoir, and other parameters are defined before.
Here, the defined optimization model is developed and extended to define an optimization model for reservoir rule curves and policies determination considering the current conditions. Generally, in this problem, the water releases values can be defined as a linear or nonlinear polynomials function of water inflow or water storage values at the different operation times and other effective parameters. In this research, two linear functions (Eqs. 10 and 11) and one nonlinear function (Eq. 12) are defined to determine the operation rule curves of the reservoir. It is worth noting that based on the previous researches, the typical patterns of the equations are defined as follows (Mousavi et al. 2007;Bozorg Haddad et al. 2008 where R y;m = water releases from the reservoir at mth month and yth year, I y;m = water inflow into the reservoir at mth month and yth year, S y;m = water storage volume of the reservoir at mth month and yth year, P y;m = the precipitation values at mth month and yth year, and N y = the number of total years, and other parameters are defined before. These newly defined Eqs. (11, 12, and 13) are replaced with water release parameters (RðtÞ) of Eq. 7, leading to three different optimization models for reservoir rule curves determination. These optimization models can be solved using different methods in which the NLP method (GAMS software) is used here to solve these optimization models.
Generally, different criteria such as root mean squared error (RMSE), correlation coefficient (R), and Nash Sutcliffe (NS) are usually used to evaluate the performance and to select the best model for water inflow prediction. These criteria are defined as (Gong et al. 2016): where X m;i = ithe computational data, X o;i = ithe observational data, X m;i = average value of computational data, X o;i = average value of observational data, and N is the number of data. It is worth noting that several indexes such as the reliability, resiliency, vulnerability, maximum deficit, and sustainability indexes have been proposed and used to evaluate the performance of different proposed models. These indexes can be defined as (Safavi et al. 2015; Sandoval-Solis et al. 2011): Res ¼ 100 1 where Rel time = reliability index of time, Rel volume = reliability index of volume, f = the total number of failures (the number of operation periods that the goal is not obtained), n = the total number of operating periods, Max Def = maximum deficit for all operating periods, Res = resiliency index, f s = the total number of continues failures, Vul = vulnerability index, SI = the sustainability index, and other parameters were defined before.

Results and discussion
In this section, the operation rule curves and policies of the Zayandehroud dam reservoir are presented and compared using different indexes such as SI. These results are determined here using SOP, SDP, NLP (solving three defined models), and both cases and approaches of ANN and GP models. It is worth noting that water loses values from the Zayandehroud dam reservoir at operation periods are assumed here to zero due to lack of data and information.
At first, the reservoir operation policies are determined using the SOP method in which these results are presented in Fig. 5. Then, the operation policies are determined using the equal number of data approach for the SDP method. Table 1 shows the obtained results when the number of intervals is equivalent to five. In this table, k denotes the reservoir storage volume interval and i denotes the water inflow interval. In addition, the infeasible solutions are presented with small horizontal dash lines. These values are corrected using the method proposed in Sect. 4. The probabilities values of Markov matrices are presented in Tables 2 and 3 considering five and ten intervals for water inflow, respectively. In addition, Tables 4 and 5 show the boundaries of the monthly state variable (water storage volumes) classes for five and ten time interval cases, respectively. Furthermore, the obtained water releases from the reservoir (operation policies) are presented in Fig. 6 using the proposed process of SDP method for five and ten time intervals cases.
Comparison of the results of the SDP method shows that the results of five time interval case are better than the ten time interval case due to the fact that the number of infeasible solutions is smaller in this case. This fact is highlighted at the end of this section by computing the water resource indexes such as SI.
In this research, the optimal operation policies of the Zayandehroud dam reservoir over 38 years are also determined by solving three proposed mathematical optimization models using GAMS software. The constant coefficients of operation policies of the first model are presented in Table 6 for all months. In addition, the constant coefficients of operation policies of the second and third models are presented in Tables 7 and 8 4 110  242  175  15  -------24   5 274  204  259  71  36  ------163  3 1 -76  88  125  136  52  85  ----114   2 71  132  110  137  145  63  96  33  -13  -139   3 294  179  142  13  ------138  reservoir (operation policies) are presented in Fig. 7 using all proposed optimization models. Comparison of the results shows that the results of the third proposed model are better than other proposed models due to smaller objective function values. In addition, the first model is a simple and applicable model with a smaller SI index value. These facts are highlighted at the end of this section by computing the water resource indexes such as SI and objective function values.
To show the efficiency of the GP model, here, this model is also used to determine the operation policies of the Zayandehroud dam reservoir by proposing two cases and approaches. For this purpose, the results of the third mathematical optimization model are used for modeling. At first, it is necessary to determine the proper values of the GP parameters in which they are obtained by the sensitivity analysis approach. In other words, these values are obtained using usual and formal methods, means trial and error method. The obtained proper values are shown in Table 9. At first, the obtained results of the first proposed cases are presented for this case. The influential input variables for each month are shown in Table 10. It is worth    noting that, in this case, the obtained objective function value and SI indexes are 33.94 and 51.3%, respectively. Comparison of the results shows that the acceptable values can be determined for water storages volumes at the end of the month using both proposed approaches of the GP model. Figure 8 shows the predicted values of water storage volumes at the end of the month using the GP model compared to actual observed values. It is worth noting that different equations can be proposed to predict the water storage volume of the Zayandehroud dam reservoir at the end of the month using the first case of the GP model. The obtained equations for each month, which can be used as reservoir operation rule curves, are presented as follows: 22 May to 21 June: S(t + 1) ¼ sinðsinðsinðsinðsinðSðtÞÞÞ^((0.36048 + I(t -5))Î(t)))))  Fig. 6 Water release values of the Zayandehroud dam reservoir using SDP (considering five and ten time intervals) where all the parameters were defined before. These equations can be used to determine the water storage values at the end of operation period (monthly operation period) based on the water inflow, precipitation, demand, and storage volume values. The operation policies of the Zayandehroud dam are determined using both approaches of the second proposed case of the GP model. The obtained results of the first approach of the second GP case are presented in Table 11, considering two conditions (with and without considering time index, T). Figure 9 shows the water storage volume values at the end of the operation period (month) using the first proposed approach of the second GP case with and without considering time index T. Furthermore, the obtained equations for reservoir operation policies determinations are presented as Eqs. 35 and 36 with and without considering time index T.
Sðt þ 1Þ ¼ sinðsinðsinðsinðsinðsinðSðtÞ þ 0:62626 Â sinðsinðIðtÞÞÞ Â sqrtðsinðsinðIðtÞÞÞÞÞÞÞÞÞÞ where all the parameters are defined before. These equations can be used to determine the water storage values at the end of the operation period (monthly operation period) based on the water inflow, precipitation, demand, and storage volume values. Comparison of the results shows that using this proposed approach leads to acceptable results for the Zayandehroud dam operation policies in which the accuracy of the results is improved by considering time index T. Finally, the obtained results of the second approach of the second GP case are presented in Table 12, considering two conditions (with and without considering time index, T). Comparison of the results shows that using this proposed approach leads to acceptable results for Zayandehroud dam reservoir operation policies in which the obtained equations are presented as Eqs. 37 and 38, respectively, with and without considering time index T. Figure 10 shows the water storage volume values at the end of the operation period (month) using the second proposed approach of the second GP case with and without considering time index T. Comparison of the results shows the better performance of first proposed approach for reservoir operation policy determination.
Sðt þ 1Þ ¼ sinðsinðsinðIðtÞ Â sinðsinðsinðSðt À 10Þ Â IðtÞ þ SðtÞÞÞÞ þ sinðsinðsinðsinðsinðsinðSðtÞÞÞÞÞÞÞÞÞÞ where all the parameters are defined before. These equations can be used to determine the water storage values at the end of the operation period (monthly operation period) based on the water inflow, precipitation, demand, and storage volume values. It is worth noting that the obtained objective function value and SI index are 33.94 and 51.3%, respectively, using the first case of the GP model. In addition, the obtained objective function value and SI index are 32 and 49.3%, respectively, using the second case of the GP model. Furthermore, Fig. 11 shows the best water storage volume values at the end of the operation period (month) using both proposed cases of the GP model.
Finally, the ANN model is also used here to evaluate the performance of proposed approaches and cases of the GP method for determining the Zayandehroud dam reservoir operation policies. For this purpose, the results of all mathematical optimization models are also used for modeling. Table 13 shows the structure of the ANN model for is worth noting that the structure of the second and third models is similar to that of the proposed structure presented in Table 13, and only the output vectors are different. The monthly water release values of the Zayandehroud dam reservoir are shown in Fig. 12 using the third proposed model of ANN. Fig. 7 Water release values of the Zayandehroud dam reservoir using all proposed mathematical optimization models  T, S(t), S(t-1), S(t-10), I(t), I(t-1) 32 49.3 S(t), S(t-1), S(t-9), I(t), I(t-1), I(t-7), I(t-9), I(t-11), P(t-1), P(t-3), P(t-4), P(t-9), D(t), D(t-1), D(t-4), D(t-5), D(t-6) 32.1 49.7 Finally, the obtained values of all defined indexes for all proposed and used methods are presented in Table 14. The obtained results indicate that the results of proposed GPbased methods are acceptable in comparison with other proposed methods. In addition, the proposed GP-based models can overcome the limitations of other proposed methods leading to simple and more applicable operation policies for the Zayandehroud dam reservoir. At the end of this research, a critical issue should be noted. Here, the reservoir operation rule and policies of the Zayandehroud dam reservoir are determined to satisfy the downstream water demand without considering social issues such as farmers' behavior in anthropogenic droughts condition. Therefore, extending and developing the proposed methods to provide future insights for policy-makers and assess the performance of complex water resources systems such as the methods of Pouladi et al. (2019 and is an essential issue, and therefore, it should be considered for future studies.

Conclusion remarks
In this research, two different approaches for two suggested cases were proposed to determine the operating policies of the Zayandehroud dam reservoir using GP. In the first approach, by hybridizing the ANN model with the NSGA-II algorithm, a hybrid method was used to reduce and select influential input data set. However, in the second approach, the influential input variables were chosen automatically using the GP method. The obtained results were compared with the SOP, SDP, ANN model, and NLP proposed methods by computing the reliability, resiliency, vulnerability, maximum deficit, and sustainability indexes. For NLP methods, here, three different mathematical optimization models were proposed and solved using GAMS software. Comparison of the results showed that using the SDP method led to an infeasible solution, which means inaccurate value for reservoir operation policies of some months. However, this problem was solved by considering the values of allowable water storage volume bounds and computing new values for the water release. In addition, for the best ANN model, the RMSE and R 2 values of the training (test) data set were 0.33 (10.6) MCM and 1 (0.98), respectively, using results of the third mathematical optimization method. Furthermore, by using the results of the third mathematical optimization method, the second proposed GP case led to simple and more applicable reservoir operation policies in comparison with the first proposed case. In other words, the obtained water deficit values of the second proposed GP case were smaller than the first one with the value of 32, in which in this case, the SI value was 49.3%. Finally, a comparison of the obtained indexes values of all proposed and used methods showed that the results of proposed GP-based methods are acceptable in comparison with other proposed methods leading to simple and more applicable operation policies for the Zayandehroud dam reservoir. For future studies, extending these methods based on the socio-hydrological modeling framework for assessing the performance of complex water resources systems should be noticed.
Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sector.

Declarations
Conflict of interest Both authors declare that he and she have no conflict of interest.
Data availability My manuscript has associated data in a data repository.
Ethical approval This article does not contain any studies with human participants or animals performed by any of the authors.