Comparative Study of Optimization Algorithms for The Optimal Reservoir Operation

Apart from traditional optimization techniques, e.g. progressive optimality algorithm (POA), modern intelligence algorithms, like genetic algorithms, differential evolution have been widely used to solve optimization problems. This paper deals with comparative analysis of POA, GA and DE and their applications in a reservoir operation problem. The results show that both GA and DES are feasible to reservoir operation optimization, but they display different features. GA and DE have many parameters and are difficult in determination of these parameter values. For simple problems with mall number of decision variables, GA and DE are better than POA when adopting appropriate parameter values and constraint handling methods. But for complex problem with large number of variables, POA combined with simplex method are much superior to GA and DE in time-assuming and quality of optimal solutions. This study helps to select proper optimization algorithms and parameter values in reservoir operation.


INTRODUCTION
Reservoir management and operation are one of the most complex problems in water resources management (Simonovic 1987), primarily due to (1) the stochastic nature of stream flows, (2) the multi-objective nature of reservoir operation, (3) the multistage and dynamic nature of reservoir operation decision-making. It is very difficult to find a single model or technique that is universally accepted for operation of reservoir systems (Liu et al. 2011). Due to complexity of reservoir operation problems, lots of optimization algorithms have been introduced and applied to reservoir management and operation, including some traditional methods and modern intelligence algorithms. The performance of these algorithms also carried out a number of comparative studies based on some standard test functions or highly simplified engineering problems; however, these testing problems cannot reflect the complexity and diversity of the calculation of the actual engineering applications. For the optimal reservoir operation, there is still a lack of detailed and comprehensive discussion on the performance of algorithms. In this paper, some optimization algorithms, including traditional algorithm (e.g. POA) and modern intelligent algorithms (e.g. GA and DE), are selected to have a detailed comparative analysis on the performance of reservoir operation from aspects of parameter, operator, variable scale, and constraint handling methods, etc. The results can help to the selection of optimization algorithms for the reservoir operation.

DETERMINISTIC OPTIMIZATION MODELS OF RESERVOIR OPERATION AND CONSTRAINT HANDLING METHODS
As an illustrative example, China's Three Gorges Reservoir (TGR) is selected as a case study to evaluate the performance of optimization algorithms. The TGR is a quarterly regulated multipurpose reservoir, with a length of 660 km and a gigantic flood storage capacity of 22.15 × 109 m 3 , and plays a very important role in the flood control of the Yangtze River. The practical operation of the TGR needs to consider many utilization objectives, such as flood control, power generation, shipping transport and so on.

Deterministic Optimization Models of Reservoir Operation.
To facilitate comparison of the algorithm performance, flood control is selected as the operation objective. For an inflow flood, the operation goal of the TGR is to discharge flood water as much as possible on the premise of ensuring the flood control safety of the dam and downstream, in order to vacate enough capacity for the following floods. In the flood control operation model, the decision variable is outflow at each operation stage, and the objective function is to maximize the average outflow with the constraints of mass balance requirement, reservoir storage limits, power generation limits and release requirement, i.e. 1 1 max where T is the total number of time periods in calculation, and t Δ is the length of time for each period.
A is the coefficient of hydropower generation. t The 100-year design flood based on the 1954 typical flow is selected to evaluate the performance of optimization algorithms. The design flood lasted 30 days, and covered 120 times periods with 6 hours an calculation period. With n defined as the number of continuous periods from the first period, several optimization problems with different numbers of decision variables can be designed when n takes different values, such as, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110 and 120.

Constraint Handling Methods.
Complex constraints handling has been the emphasis and difficulty in the practical application of the various algorithms, and have a direct influence on the performance of algorithms. There are many constraints handling methods, such as discarding infeasible solutions, repairing infeasible solutions, penalty function, multiobjective method, and so on (Coello 2002), among which repairing infeasible solutions and penalty function are more feasible methods for the high-dimensional multiconstraint problems of reservoir operation. Repairing algorithm is problem-dependent, and needs to be designed according specific reservoir operation problem. Penalty function approach has to be introduced extra empirical parameters to the violation penalty of infeasible solutions to find the intersection of a line with the boundary of the feasible region. Howson and Sancho (1975) presented the progress optimality algorithm (POA) for the solution of multi-state dynamic programming problems. In the algorithm, progressive optimality theory based on Bellman's principle of optimality was introduced, which indicates that each pair of decision sets is optimal in relation to its initial and terminal values (Howson 1973). Traditional algorithm is conceptually simple, iteratively applying a general two-stage solution to successive overlapping stages of the problem, and giving successive approximations that converge to the optimal path (Howson 1975). We consider that the objective function is

PROGRESS OPTIMALITY ALGORITHM
where X 0 , X N are given. This is the equivalent to the twopoint boundary value problems in dynamic programming. The algorithm solves the minimization of a two-stage problem where Xj, j=1,2,3...N are vectors of stage j (Howson 1975).
At present research, the ways of discrete domain have been widely used to get an optimal value. However, it limits the accuracy and capacity of optimal search in some way. So this paper adopts simplex method to deal with two-stage problems of decision optimization search. Simplex method has been presented originally by Splendy et al. (1962). And then Nelder and Mead (1965) have done some improvements aiming at difficulties that they have met in accelerating search or searching in the valley and ridge. The improved method allow simplex to vary its shape, size and orientation to adapt itself to the local contour of the objective function (Shu et al. 1957). The method described for the minimization of a function of n variables, which depends on the comparison of function values at the (n+1) vertices of a general simplex, followed by the replacement of the vertex with the highest value by another point. The simplex adapts itself to the local landscape, and contracts on to the final minimum (Splendy et al. 1962).

GENETIC ALGORITHM
Genetic Algorithm (GA) was first introduced by Holland in 1975 .GA is a search heuristic to generate solutions to optimization problems using techniques inspired by natural evolution, such as inheritance, mutation, selection, and crossover. The main characteristics of GAs include: 1) In the genetic algorithm, it is a practical matter about the components of genetic chromosome which need encoding firstly and replace the original mathematical expression (Yan et al. 2016). Therefore, GAs can solve discrete, non-convex, discontinuous problems without differentiation. 2) These methods are inherently parallel, The simple GA, mainly adopting binary coding, simple arithmetic crossover and uniform mutation, has significant limitations in practical engineering problems. The crossover operator is the main operator and has a strong global search capability. For real-coded genetic algorithms, there have been proposed many crossover methods by improving simple arithmetic crossover, including: blend crossover (BLX) (Takahashi and Kita 2001), simulated blend crossover (SBX) (Deb and Agrawal 1995), unimodal normal distribution crossover (UNDX) (Ono et al. 1999) and simplex crossover (SPX) (Tsutsui 1999). The blend crossover (BLX) combines parental information in a component-wise manner in generating offspring and shows good search ability in optimizing separable objective functions (Wen 2006). The mutation operator is an auxiliary operator in genetic algorithm and has certain local search capability. Several mutation operators are widely used in real-coded genetic algorithm (e.g., uniform mutation operators, boundary mutation operators, non-uniform mutation operator, Gaussian mutation operators). The non-uniform mutation operator  is relatively simple and efficient, and its search space gradually narrowed with the increase of generation. This paper mainly adopts a more efficient GA, which uses blend crossover (BLX-α) and non-uniform mutation, and analyzes its performance in optimal reservoir operation. α is a constant on the interval (0,1). In addition, the non-uniform mutation contains a parameter β with the integer value ranging from 2 to 5, which determines the degree of dependence of random number perturbation on generation. Both α and β need to be selected based on experience.

DIFFERENTIAL EVOLUTION
Differential Evolution (DE) is a very simple population based, stochastic function minimizer which is very powerful at the same time. It was first appeared as a technical report in 1995. DE operates through similar computational steps as a standard GA. However, unlike traditional GA, the defining characteristic of DE is the generation of offsprings (the so called trial vector) by employing the differential mutation operation followed by crossover of mutant vector with a predetermined parent vector (Thangavelu and Velayutham 2015). DE generates new parameter vectors by adding the weighted difference between two population vectors to a third vector. If the trial vector yields a lower cost function value than the target vector, the trial vector replaces the target vector in the following generation. The detail of DE is described by Price et al. (2006).
In order to classify the different kinds of DE model, the notation DE/x/y/z, is introduced here. x specifies the vector to be mutated which currently can be "rand" (a randomly chosen population vector) or "best" (the vector of lowest cost from the current population); y is the number of difference vectors used; z denotes the crossover scheme. The current variant is "bin" (Crossover due to independent binomial experiments) .
There are two main control parameters of DE: scale factor F and crossover rate CR. F is a constant number ranging from 0 to 2, which controls the diversity and convergence of the population by controlling the scaling degree of difference vector. And CR is a user chosen parameter on the interval (0, 1), which controls the degree of participation of each dimension in the individual parameters and the balance between global and local search ability.

GA AND ITS APPLICATION TO OPTIMAL RESEVOIR OPERATION
GA is a search heuristic to generate solutions to optimization problems using techniques inspired by natural evolution, such as inheritance, mutation, selection, and crossover. GA was described by . The simple GA, mainly adopting binary coding, simple arithmetic crossover and uniform mutation, has significant limitations in practical engineering problems. In this paper, a more efficient GA is introduced, which uses blend crossover (BLX) and non-uniform mutation based on floating-point encoding. The detail of BLX and non-uniform crossover are described by  and Takahashi and Kita (2001), respectively.
There are several important parameters for the running of the proposed GA. The parameters α in the BLX and β in the non-uniform mutation are set to be 0.7 and 2, respectively, which was found to be more efficient by the author. For the problems with different number of variables, different crossover rates (Pc) and mutation rates (Pm) were tested in this study. The population size is 500, and the optimization process terminates when the optimal function value achieves the global optimal solution, the maximum generations is more than 50000, or the optimal function value changes less than 1.0. For each parameter pair of the decision variable number, the crossover rate and the mutation rate, the optimization program will run 10 times, of which the average value of the optimal function values and the generations are calculated as the optimal results of each parameter pair. All optimization programs are carried out on a computer with the CPU of Intel Core 2 Quad Q8200 Processor (2.33GHz) and the memory of 4GB.
The results of the problems with three different numbers of decision variables are shown in Fig. 1, Fig. 2 and Fig. 3, respectively. Fig.1 (a), Fig.2 (a) and Fig.2 (a) shows the best objective function value derived by GA with different crossover and mutation rates when the number of variables are 10, 30 and 50, respectively. Fig.1  (b), Fig.2 (b) and Fig.2 (b) shows the maximum generations with different crossover and mutation rates when the number of decision variables is 10, 30 and 50, respectively. In these figures, the parameter pair of Pc and Pm with darker color block is better. The results show that the optimal function value and the maximum generation are more sensitive to Pm. Pc has more influence on the maximum generation but little to the optimal function value. In addition, the more decision variables, the smaller the probability of achieving global optimal solution. Larger mutation rate may need more generations but achieve worse solutions. It is generally recommended to not more than 0.1 for the mutation rate, and it should be smaller when the decision variables are more. The crossover rate is preferably from 0.6 to 0.8 for 30 decision variables, and smaller for 50 decision variables.

DE AND ITS APPLICATION TO OPTIMAL RESERVOIR OPERATION
DE, developed by , is also a stochastic population-based optimization algorithm, simple yet powerful. DE generates new parameter vectors by adding a weighted difference vector between two population members to a third member. If the resulting vector yields a lower objective function value than a predetermined population member, the newly generated vector will replace the vector with which it was compared in the following generation. The comparison vector can but need not be part of the generation process mentioned above. In addition the best parameter vector is evaluated for every generation in order to keep track of the progress that is made during the minimization process. The detail of DE is described by Price et al (2005).
There are several variants of DE, among which DE/BEST/1/EXP is more efficient for the reservoir operation problem in this paper. The differential weighting factor (F) and the crossover probability (CR) are two main parameters. F is a positive real number that controls the rate at which the population evolves. While there is no upper limit on F, effective values are seldom greater than 2.0. CR is a user-defined value that controls the fraction of parameter values that are copied from the mutant with a value of 0 to 1. Other parameter values are set to be same to the GA above. For the problems with three different numbers of decision variables, the optimal global solutions are all achieved by DE. Fig.4 shows the average generation for each parameter pair. The parameter pair with darker block needs less generation and so is better. The results show that the performance of DE with small value of F and large value of CR is better. However, the performance may be poor when CR is more than 0.9. With the increase in the number of decision variables, F need to be increased and CR decreased gradually.

COMPARATIVE ANALYSIS OF OPTIMIZATION ALGORITHMS FOR THE RESERVOIR
GA and DE can achieve best results in application to reservoir operation, yet needing proper parameter values. Here, the performances of these modern intelligent algorithms are compared with a traditional algorithm POA.

POA in combination with simplex method.
POA is presented by Howson and Sancho (1975) for the solution of multi-state dynamic programming problems. It is a method of successive approximation using a general two-stage solution. The algorithm is computationally efficient and has minimal storage requirements. For each two-stage solution, simplex method for function minimization is used. Simplex method was introduced by Splendy et al (1962) for tracking optimum operating conditions by evaluating the output from a system at a set of points forming a simplex in the factor-space, and continually forming new simplexes by reflecting one point in the hyperplane of the remaining points. The method was improved by Nelder and Mead (1965). The improved simplex adapts itself to the local landscape, and contracts on to the final minimum. The method is shown to be effective and computationally compact.

8.2
Parameter setting of optimization algorithms.
POA adopts penalty function in constraint handling. GA and DE adopt penalty function and repairing algorithm for unfeasible solutions. Parameters of GA and DE are all determined with quite good performance. The crossover and mutation rate are 0.70 and 0.05, respectively. The differential weighting factor and the crossover probability based on DE/BEST/1/EXP are 0.4 and 0.8, respectively. The population sizes of GA and DE are all 500, and the optimization process terminates when the optimal function value achieves the global optimal solution, the maximum generations is more than 50000, or the optimal function value changes less than 1.0. For each number of the decision variables, the optimization program will run 10 times, of which the average value of the optimal function values and the generations are calculated as the optimal results. Generation of POA is defined as the number of loop iterations.

8.3
Comparative analysis of optimization algorithms.
The results are shown in Table 1. It can be seen that POA combined with simplex method is superior to GA and DE combined with penalty function in performance. For each number of decision variables, POA can achieve the global optimal solutions. The performance degrades little with the increase of decision variables, while the probability to the global optimal solutions of GA and DE decreased obviously and time-assuming increases significantly with the increase of variables. Compared with GA, DE combined with penalty function has better solutions and less time-consuming, but hard to obtain feasible solutions for more than 70 decision variables. In combination with GA, repairing algorithm can decrease time-assuming greatly when the number of variables is smaller than 90, but may achieve poorer solutions when then number of variables is larger than 40. In combination with DE, repairing algorithm can also decrease time-assuming greatly and achieve better solutions than penalty function. Repairing algorithm is more efficient than penalty for DE without the limitation of number of variables, but for GA only when the number of variables is small. In addition, GA and DE combined with repairing algorithm are better than POA when the number of variables is less than 30.

CONCLUSIONS
GA and DE are all feasible in application to reservoir operation, but have many parameters and are difficult in determination of these parameter values. For simple problems with small number of decision variables, GA and DE are better than traditional algorithms (e.g. POA) when adopting appropriate parameter values and constraint handling methods. But for complex problem with large number of variables, POA combined with simplex method are much superior to GA and DE in timeassuming and quality of optimal solutions.