Comparison of representative heuristic algorithms for integrated reservoir optimal operation

Heuristic algorithms (HAs) are widely used in integrated reservoir optimal operation due to the fast calculation and simple design. Existing literature usually focuses on one or two categories of HAs and simply reviews the state of the art. To provide an overall understanding and specific comparison of HAs in integrated reservoir optimal operation, differential evolution (DE), particle swarm optimization (PSO), and artificial physics optimization (APO), which serves as typical examples of the three categories of HAs, were compared in terms of the development and applications using a designed experiment. Besides, the general model with constraints and fitness function, and the solving process using a hybrid feasible domain restoration method and penalty function method were also presented. Taking a designed experiment with multiple scenarios, four evaluation criteria are used for the comprehensive comparison. Results of the comparison show that (a) the problem of integrated reservoir optimal operation is a mathematic programming problem with the narrow feasible region and monotonic objective function; (b) it is easy to obtain the same optimal objective function value but different optimal solutions when using heuristic algorithms; and (c) DE, PSO and APO does not result in a clear winner, but DE could be more appropriate for integrated reservoir operation optimization according to the four evaluation criteria used in this study.


Introduction
The reservoir is the most effective engineering infrastructure that provides water storage and redistributes water inflow in time and space. Most reservoirs, especially large-scale water management projects, are designed as integrated projects that may provide for flood control, municipal and industrial water supply, hydroelectric power, agricultural irrigation, navigation, and wetland protection (Labadie 2004). Generally, reservoir operation is trying to answer the following question: when to release the water from the storage and how much water should be released at that time? This is a typical sequential decision-making problem (Oliveira and Loucks 1997). However, various factors and restrictions involved in the problem, such as politics, emergency, precipitation monitoring, flood forecasting, state of the equipment, water demand, decision-makers' preferences, uncertainty make the real-time reservoir operation quite complicated (Rani and Moreira 2010). Optimization techniques of integrated reservoir operation remain one of the hot topics in water resources management.
Optimization techniques have been gradually introduced to reservoir management and operations since the 1960s. Then a large number of publications and achievements have emerged with the classical optimization techniques, which is analytical. Yeh (1985) reviewed the reservoir management and operations models, including linear programming (LP), dynamic programming (DP), and non-linear programming (NLP). Generally, the traditional techniques contribute to the optimization with clear mathematical formulation, objective function, and constraints, but also exhibit limitations on non-convex optimization problems and real-time reservoir operation implementation (Ateş and Yeroğlu 2018). Besides, the computational cost of classical techniques grows considerably as the reservoir number increases (Sharif and Swamy 2014).
As computational intelligence rapidly develops in the past two decades, heuristic algorithms based on computational intelligence are becoming significant approaches for solving optimization problems and widely explored in solving reservoir operation problems (Al-Jawad and Tanyimboh 2017; Jia et al. 2015;Yang et al. 2020).
These algorithms provide fast computation, simple design, global optimal solutions, and superior generality. They can easily overcome lots of limitations of classical reservoir optimization approaches, such as equality constraints, high-dimensional computation, nonlinear relationships, calculation efficiency, multiple objectives (Dahmani and Yebdri 2020).
There are no free lunch theorems for optimization (Wolpert and Macready 1997).
Every heuristic algorithm has advantages but also disadvantages and different heuristic algorithms exhibit different performance on different problem instances (Chen et al. 1999). Much literature pays effort into the review of the applications of different algorithms. Labadie (2004) assessed the state-of-the-art in the optimization of reservoir system management using evolutionary algorithms, along with the application of neural networks and fuzzy rule-based systems for inferring reservoir system operating rules. Hossain and El-shafie (2013) discuss the performance of the evolutionary algorithms and the ability to integrate with other techniques on the ground of the benchmark problem. The results verify that the evolutionary algorithms overcome the drawbacks of traditional algorithms and the swarm intelligence is less comparatively than GA but provides a great scope for further development. Maier et al. (2015) reviewed the potential applications of different evolutionary algorithms in various application areas and paid significant research effort on the better solutions of reduced computational effort.
However, existing reviews usually focus on one or two types of heuristic algorithms comparing the different variants, which can not give an overall understanding of the advantages and disadvantages among the three categories. Besides, little literature gives a particular comparison of their performance of multi-objective reservoir optimal operation. To bridge these gaps, this paper presented specific comparison research with a designed experiment of differential evolution (DE), particle swarm optimization (PSO), and artificial physics optimization (APO) algorithms, which respectively represent evolutionary computation algorithms, swarm intelligence algorithms, and physical and chemical law-based algorithms. The goals include (a) characterizing multi-objective reservoir management and operation optimization problem; (b) identifying the difference between three typical heuristic algorithms on multi-objective reservoir optimal operation; (c) explaining why the three selected heuristic algorithms result in different performance.

Differential evolution
Differential evolution (DE) was originally designed to solve the Chebyshev Polynomials and then widely used in reservoir operation. The core components of the DE algorithm include mutation operator, crossover operator, and selection operator.
(1) Mutation operator DE can be divided into different modes based on the numbers and types of the selected parent population. The general convention used is DE/α/β/γ, where α is the vector to be perturbed (can be Rand: randomly selected or Best: best population), β the number of vectors considered for the mutation and γ the type of crossover used (can be Exp: exponential or Bin: binomial). The most commonly used two types of mutation operator are as follows: where Xt () is an individual; PS is the parameter scaling factor that controls the scaling weighting in mutation operator; i is an identifier for one individual; P N is the total number of individuals in a population; g is an identifier for evolutionary generations; best is an identifier for the individual that has the best fitness; r1, r2, and r3 are random integers located in the interval [1, P N ], they are not equal to each other as well as i , and also they are used as identifiers for individuals.
(2) Crossover operator Test population will arise from the use of the crossover operator. The test population is generated between the parent population and mutation population according to the following (for minimization optimization problem): Yt () is a test individual; PC is the crossover rate parameter that controls the probability of mutation-to-test individual in crossover operator; Rand is a random value within uniform distribution [0, 1]; other variables as shown above.
(3) Selection operator DE algorithm is a "greedy" selection, namely: comparison the fitness function value of test population with parent population, then the outstanding population is preserved and into the next evolution of the offspring population. The selection operator can be expressed as for minimization optimization problem): where  F() is the fitness function which is used for assessing the merit of the individual; other variables as shown above.
Since DE was proposed, many applications were documented in the literature as Table 1 shows. DE is computationally very efficient and offers a high level of robustness. Besides, DE is one of the best evolutionary algorithms for solving optimization problems with real-coded variables, which makes it more attractive for implementation in practice. Table 1. Review of reservoir operation using differential evolution Paper Research Reddy and Kumar (2008) adopted multi-objective DE for the simultaneous evolution of optimal cropping pattern and operation policies for a multi-crop irrigation reservoir system with the objectives that maximize total net benefits and the total irrigated area. Regulwar et al. (2010) presented DE for the optimal operation of the multi-objective reservoir, in which the objective was the maximization of the hydropower production, while the irrigation supply was considered as a constraint. They showed that DE/best/1/bin was the best strategy to obtain the optimal solution.

Particle swarm optimization
Particle swarm optimization (PSO) is based on the foraging behavior of birds and takes advantage of sharing information (velocity and position) among particles (vivid birds) in the whole swarm (vivid bird group). The movement of the whole swarm is generated from disorder to order evolutionary process in a search space of the optimization problem. The optimal solution is obtained within a finite number of evolutions.
The basic components of PSO are velocity and position operator. Each individual can be seen as a particle, which has no weight and volume, while a particle flies at a certain velocity in search space. The flight velocity is dynamically adjusted by the flight experience of the individual and the flight experience of the group. The equations of velocity and position operators for updating particles are as follows: where Xt () is a position vector, which represents a particle (an individual); Vt () is a velocity vector, which represents an attribute in particle; w is the parameter of inertia weighting; c1 is the parameter of social factor; c2 is the parameter of cognitive factor; r1 and r2 are two random values within uniform distribution [0,1]; gbest is an identifier for the particle that has the best fitness in all past g populations; pbest is an identifier for the particle that has the best fitness in g-th population; other variables (i, g, NP) are as shown above in the discussion of DE.
Different variants of PSO were proposed as Table 2 shows. The PSO is characterized by the easy implementation, high accuracy, and rapid convergence.
Therefore, it attracted the attention of water resources specialists and demonstrated its superiority in solving water resources management problems. Table 2. Review of reservoir operation using particle swarm optimization

Paper Research
Clerc and Kennedy (2002) presented a constriction factor χ to replace the parameter of inertia weight w to ensure the PSO convergence. Ratnaweera et al. (2004) presented a modified PSO by the introduction of time-varying acceleration coefficients: TVAC. Kumar and Reddy (2007) adopted PSO to derive reservoir operation policies for multipurpose reservoir systems. The PSO algorithm was further improved by incorporating a new strategic mechanism called elitist-mutation to improve its performance in their study.

Mandal and
Chakraborty (2011) optimized the short-term combined economic emission scheduling proposed a hybrid multi-objective PSO-EDA algorithm, in which divided the particle population into several sub-populations and builds probability models for each process, for solving single reservoir multi-objective flood control problem.
Niu et al. (2018) came up with a parallel multi-objective particle swarm optimization algorithm, in which the swarm with large population size was divided into several smaller subswarms to be simultaneously optimized by different worker threads, for solving cascade hydropower reservoir operation balancing benefit and firm output problem.

Dahmani and
Yebdri (2020) presented a hybrid algorithm of particle swarm optimization and grey wolf optimizer for solving a single reservoir long-term operation problem with the objectives to minimize water supply deficits, results showed better performance in overcoming trapping in the local minima. (1) Mass operator

Artificial physics optimization
The mass of an individual is not a constant, but a function related to the definition of fitness value of an individual, and this mass function must be non-negative, bounded, and monotonic. The value of the mass function should be restricted between 0 and 1, and can be used as the measure of fitness value. The bigger the value of mass of the individual is, the better the fitness of the individual is. Therefore, the equation of the mass operator derived from the concave function is: where M is the mass of an individual;  F() is the fitness function which is used for assessing the quality of an individual; Xt () is a position vector, which represents an individual; exp is natural exponential; best and worst are identifiers for the individuals that have the best and worst fitness values, respectively; other variables (i, g, NP) are as presented above in the description of DE.
(2) Force operator The rule of gravitational force is: an individual with better fitness attracts the individual with poor fitness; the individual with the optimum fitness level has an absolute attraction, and attracts all the other individuals not affected by any force.
Taking the minimization problem as an example, the force operator among individuals is: where CF t () is the virtual component force vector between i-th and j-th individual; G is the gravitational constant; Rt () is the distance vector, which can be calculated as where w is the weighting parameter of inertia; r is a random variable, which obeys the uniform distribution [0,1] or the normal distribution (0,1); and other variables as shown above.
The review of APO is given in Table 3. As Table 3 shows, APO is seldom used in reservoir operation.

Integrated reservoir optimal operation model
The main purpose of the multi-objective reservoir optimal operation is finding the best way of meeting the water demand of multiple water users. The water users may include flood control, municipal and industrial water supply, hydroelectric power, agricultural irrigation, navigation, etc. When developing the long-term optimization operation model for a multipurpose reservoir, important purposes such as hydroelectric power generation and water supply are given priority (incorporated as the objective functions), and other purposes are considered as constraints.

Objective functions
For multi-objective reservoirs, with hydropower generation and water supply as the main users, maximization of hydropower generation and minimization of water supply deficit are the two most common goals considered in reservoir optimal operation.
In this paper, the supplied water refers to the water demand taken directly from the reservoir by the diversion canal which does not participate in hydropower generation.
The two goals can be expressed as the following objective functions: where B is the benefit function; C is the cost function; T is the total period of operation;  is the hydropower efficiency coefficient; t H is the hydraulic head during the period t; t HR is the power release during the period t ; t ID is the water supply-demand (such as agricultural irrigation and domestic water) during the period t ; t IR is the water supply from the reservoir during the period t.
The two objectives conflict with each other. Maximization of hydropower production requires a higher volume of water stored in the reservoir to produce more electricity, while minimization of water supply deficit tries to minimize the water supply deficit requiring more water from the reservoir to meet the demand of water supply. To handle this multi-objective competition problem, a standardization approach was applied. It is able to express two objectives in the same units by treating the water supply deficit ( where Φ is the objective function of the reservoir optimal operation model.

Constraints
In the multi-objective reservoir optimal operation model, the following constraints are considered.
(1) Water balance: t L is the total water loss (such as evaporation and leakage) during the period t; and other variables as discussed above.
(2) Water storage limits: where t V ,min and t V ,max are the lower and upper bounds, respectively, of water storage during the period t .
(3) Water for power generation limits： where H min is the minimum hydraulic head.

Reservoir optimization with heuristic algorithms
The common heuristic strategy is usually adopted when heuristic algorithms are In this paper, a general process of implementing a heuristic algorithm to solve the multi-objective reservoir optimal operation model is put forward by using a hybrid feasible domain restoration method and penalty function method. The comparison of the three selected heuristic algorithms is possible within the general process presented in the following section.

Process of computation
Step 1. Individual initialization A set of randomly generated reservoir outflows or reservoir water levels is regarded as an individual (for example, the individual will be represented by the reservoir outflow sequence in this paper), and the reservoir outflow sequence is encoded real numbers corresponding to each attribute. In this way, the search space of the heuristic algorithm is covered by randomly distributed multiple individuals that form a search swarm.
Step 2. Search space definition It is necessary to make clear the relationship between the feasible domain and the infeasible domain due to the multiple constraints in the reservoir optimal operation model. In the model, an individual of the heuristic algorithm is selected to be water release for hydropower generation ( t HR ) and water release for water supply ( t IR ) .
Therefore, the water outflow and water level can be calculated using the water balance constraint (Eq. 13) together with the initial conditions (Eq. 18). The feasible domain restoration method is adopted for the power generation limits (Eq. 15) and the water supply limits (Eq. 16). This process needs a restoration operator (see Section 4.2). The violations between the remaining unrestored constraints (Eq.14, Eq.17, and Eq.19) and individual attributes are calculated separately. The results are collectively calculated using the violation function of the other constraints ( ( )  Viol ).
Step 3. Fitness space definition A clear relationship between the objective function domain and the non-objective function domain in the multi-objective reservoir optimal operation model is necessary. The fitness space is formed when the penalty function method is implemented to establish a penalty function to be added to the objective function (Eq. 12). This process is called fitness function design (see Section 4.3). Then, following the principle of the "bigger the better" fitness function value, the individual carrying the maximum fitness function value is identified.
Step 4. Generate new individuals Various operators, such as mutation, selection, position, mass, etc., are used to update the parental individuals, and to generate new offspring individuals for the next assessment and iteration. (Different operators mat carry different names in different heuristic algorithms, however, they have the same meaning).
Step 5. Termination condition In this step, the check is performed to see whether the iteration number reached the preselected maximum number of iterations g max . If this criterion is not met, the process returns to Step 2 to continue the iterative calculations; if it is met, the iteration calculation ends with the optimal individual values and their fitness.

Constraints handling process
when the t IR encoded attribute of an individual violates the water supply limits (Eq. 16), the repair of the value of the attribute is made by making the attribute value equal to the value of the limit.
This restoration operator is expressed as:

Fitness function design in the heuristic algorithm
The two constraints are restored before, while the other constraints (Eq. 14, Eq. 17, and Eq. 19) are still in the unknown state. Fitness function design, which consists of the objective function and penalty function, will help in dealing with this unknown state. The fitness function is expressed in the heuristic algorithm as: where ug is a dynamic penalty function, and it can be calculated by is the violation function of the other constraints (Eq. 14, Eq. 17, and Eq. 19).

Case study reservoir
The multipurpose reservoir located in the south-west coastal region of India, where the rainfall is abundant with a tropical monsoon climate, is selected as a case study reservoir. The rainy season of this region usually covers June, July, August, and September. The catchment of this reservoir is elongated, hilly with steep slopes with an area of about 892 km 2 . The reservoir is designed for purposes of hydropower generation and agricultural irrigation water supply, and the total storage is 2797×10 6 m 3 , corresponding to the water level of 659.9 m. The dead storage is 145×10 6 m 3 , corresponding to the water level of 609.6 m. The annual water release from the reservoir for hydropower generation is expected to be 1912×10 6 m 3 , while agricultural irrigation demand is 850×10 6 m 3 . The normal reservoir water level is 657.9 m. The comprehensive efficiency of the hydropower station is η=7.61.

Experiment scenario
To investigate the difference among the three selected algorithms, we have designed 3 experiment scenarios of the reservoir operation problem. The algorithms are set with the same common parameters, the population number, and the maximum number of iterations. The algorithm-specific parameters are set as good as possible, to provide consistent conditions for achieving the superior performance of each algorithm.
Three particular inflow years are selected, dry year (2871×10 6 m 3 total annual volume), normal year (4203×10 6 m 3 total annual volume), and wet year (6556×10 6 m 3 total annual volume). For the three heuristic algorithms, the population number is set to be 60, 100, and 150, and the maximum number of iterations is limited to 500, 800, and 1200, respectively.
The proposed procedure of solving the operation model with the heuristic algorithm is implemented using the computer programming language (VB.NET) and database technology (SQL 2012). The data and parameters are stored in the database to facilitate unified input, and each solution of the algorithm is recorded into an independent calculation module. Each algorithm is executed 10 times independently.

Evaluation criteria
Most applications of the optimization approach use the optimal objective function value as the only criterion for evaluation. Mean of optimal objective function values, the standard deviation of optimal objective function values, mean of computational time and population diversity are chosen to evaluate the specific performance of the three representative heuristic algorithms.
(1) Mean of optimal objective function values where Ave is the mean of optimal objective function values and obtained from n repetitions of independent calculations.
(2) Standard deviation of optimal objective function values where g Div is the population diversity at g-th iteration in the calculation procedure.
This criterion is derived from biology. In this paper, each heuristic algorithm has its search individual, and the individuals form a population, the proportion of the population in the whole search space can reflect the evolution degree of the heuristic algorithm.
6 Results and discussion 6.1 Mathematical model of optimal reservoir long-term operation  representation is used in this paper, to better illustrate the multipurpose reservoir long-term optimal operation model solution obtained using heuristic algorithms. The year is divided into only two seasons (wet and dry), namely T=2, and then the optimization is performed under various control conditions. 3D schematic diagrams are shown in Figure 1 and Figure 2. (1) The dimensions of the theoretically feasible region with one dimension less than the number of decision variables (see the theoretical feasible solid line in Figure   1), due to the presence of initial water level and terminal water level equality constraints.
However, because of pre-assigned optimal solution accuracy in actual reservoir operation, the feasible domain is narrow (see the green area in Figure 1).
(2) That the theoretical objective function could be monotonic (see the blue solid line in Figure 2), and the theoretical optimum objective value could be unique (see the green point in Figure 2). However, because of the optimal solution accuracy in actual reservoir operation, the objective function domain is distributed in the tangential direction of the theoretical optimum objective function value (see the yellow solid line in Figure 2).
(3) That when heuristic algorithms are used to solve the problem (see the black dot in Figure 1), it is easy to obtain the same optimum objective function values, but different optimal solutions.  The following conclusions are made from Table 4 and Figure 3 to Figure 6.

Comparison of heuristic algorithms
(1) For all three algorithms, the optimum objective values are improving with the increase of the population size. This phenomenon is reflected by the mean criterion of the optimal objective function values (see column 5 in Table 1). However, as the population size increases, the computational time for locating the optimal solution increases (see column 7 in Table 1).
(2) The mean of optimal objective function values is not a sufficient criterion to determine the best algorithm (see column 5 in Table 1 and Figure 3), because the algorithm search is random and the parameters of the algorithm have an important effect.
However, from the deviation of optimal objective function value, DE shows better stability when compared to PSO and APO (see column 6 in Table 1).
(4) The optimal objective function values obtained by the three algorithms are similar, but the search process is very different (see Figure 5, basically the first 200 steps show a fast convergence, and after 200 steps, the search becomes local). The population diversity criterion for DE still showing some numerical value at the maximum iteration step, but it has almost zero value for PSO and APO (see column 8in Table 1). So the DE algorithm maintains a certain population diversity value unchanged, but PSO and APO show a decrease in value and fluctuation. Again, the APO performance is inferior to the other two algorithms (see Figure 6). Therefore, according to these analyses, DE shows a much higher potential of not falling into the trap of local optimal solution, while PSO and APO show the opposite.

Reservoir operation comparison
The optimal solutions of reservoir operations obtained using three heuristic algorithms are shown in Figure 7. The following observations are obtained from the analysis of results in Figure 7. The results show that: (1) The optimal solution obtained by APO tends to take provide more reservoir release for irrigation (the green triangles in Figure 7). The optimum solution obtained by PSO tends to release more water for hydropower generation (the red squares in Figure 7). The DE (blue diamonds in Figure 7) performance is superior to the other two algorithms.
(2) The PSO algorithm shows higher instability and may end up with the local non-convergence solution (in Figure 7, many red squares are out of a linear relation track).
(3) The reservoir water elevation, namely the optimum solution, is significantly different from month to month. This is also indicative of the same optimum objective value and different optimal solutions.
(4) The DE algorithm could be more appropriate for reservoir operation considering all of the four evaluation criteria.

Conclusions
This paper reviewed the application of the representative heuristic algorithm of the three categories, including DE, PSO, and APO, to reservoir optimal operation and compared the specific performance with a designed experiment. A general solution process with the implementation of heuristic algorithms is developed for the problem of integrated reservoir optimal operation together with the necessary techniques for dealing with constraints and fitness function design. Four evaluation criteria are used for the comparison of three algorithms of the designed experiment. The main conclusions obtained from the comparative analyses are as follows.
(1) The problem of integrated reservoir optimal operation is a mathematical programming problem characterized by the narrow feasible region and monotonic objective function. When this problem is solved using the heuristic algorithm based on the random search method, it is easy to obtain the same optimum objective function values but very different optimal solutions. reservoir operations are taken into consideration, the DE seems to be more appropriate for integrated reservoir optimal operation than PSO and APO.