Optimization of Water-Supply and Hydropower Reservoir Operation Using the Charged System Search Algorithm

The Charged System Search (CSS) metaheuristic algorithm is introduced to the field of water resources management and applied to derive water-supply and hydro-power operating policies for a large-scale real-world reservoir system. The optimum algorithm parameters for each reservoir operation problems are also obtained via a tuning procedure. The CSS algorithm is a metaheuristic optimization method inspired by the governing laws of electrostatics in physics and motion from the Newtonian mechanics. In this study, the CSS algorithm’s performance has been tested with benchmark problems, consisting of highly non-linear constrained and/or unconstrained real-valued mathematical models, such as the Ackley’s function and Fletcher–Powell function. The CSS algorithm is then used to optimally solve the water-supply and hydropower operation of “Dez” reservoir in southern Iran over three different operation periods of 60, 240, and 480 months, and the results are presented and compared with those obtained by other available optimization approaches including Genetic Algorithm (GA), Ant Colony Optimization (ACO), Particle Swarm Optimization (PSO) and Constrained Big Bang–Big Crunch (CBB–BC) algorithm, as well as those obtained by gradient-based Non-Linear Programming (NLP) approach. The results demonstrate the robustness and superiority of the CSS algorithm in solving long term reservoir operation problems, compared to alternative methods. The CSS algorithm is used for the first time in the field of water resources management, and proves to be a robust, accurate, and fast convergent method in handling complex problems in this filed. The application of this approach in other water management problems such as multi-reservoir operation and conjunctive surface/ground water resources management remains to be studied.


Introduction
Optimal utilization of available fresh-water resources requires development and application of robust and effective methods to plan and operate current reservoirs [1,2].Effective reservoir operation requires policies that optimize releases from the reservoir or storage volume, in order to achieve desired objectives such as maximizing power generation or minimizing water deficit, flood risk, and operation costs.Adoption of these policies are a must for the operators to make decisions based on the current and past conditions of the reservoir storage and river inflow, in order to manage the upcoming flood risks and water shortages [1][2][3].
Effective management of current and planned reservoirs will become even more challenging as the Earth's climate has changed in the past and is projected to change even more in the future decades,

Charged System Search (CSS) Algorithm
The CSS algorithm is developed based on the governing laws of electrostatics in physics and motion from the Newtonian mechanics [41].The CSS utilizes a number of agents or solution candidates, which are called charged particles (CPs).Each CP is considered to be a charged sphere which can impose electrical forces on the other CPs according to the Coulomb and Gauss laws of electrostatics.Then, the Newton's law is utilized to calculate the acceleration value based on the resultant force acting on each CP.Finally, utilizing the Newtonian mechanics, the position of each CP is determined at any time step based on its previous position, velocity and acceleration in the search space [41].Each CP is considered as a charged sphere with radius a, which has a uniform volume charge density (q i ) equal to: where fitbest and fitworst are the best and the worst fitness values of all the particles, and fit(i) is the fitness of the particle i, and N is the total number of CPs.The initial positions of CPs are assigned randomly in the search space, within the boundaries determined by the problem.The initial velocities of the CPs are taken as zero.
The CPs are scattered in the search space and can impose electric forces on the others.The magnitude of the force for the CP located inside or outside of the sphere are determined differently.The resultant electrical force acting on CPs inside or outside of the sphere is determined using: where F j is the resultant force acting on the jth CP. r ij is the separation distance between two particles defined as: where X i and X j are the positions of the ith and jth CPs, respectively; X best is the position of the best current CP.ε here is a small positive number to avoid singularity.The p ij determines the probability of moving each CP toward the others as: As shown in Equation ( 2), the force imposed on a CP inside the sphere is proportional to the separation distance between particles.However, for the CPs located outside the sphere, it is inversely proportional to the square of the separation distance.The new locations of the CPs are calculated based on the resultant forces and the laws of the motion.At this step, each CP moves towards its new position according to the resultant forces and its previous velocity as: V j,new = X j,new − X j,old ∆t where rand j1 and rand j2 are two random numbers uniformly distributed in the range (0,1).Here, m j is the mass of the j th CP, which is set to be equal to q j .∆t is the time step and is set to unity.k a is the acceleration coefficient; k v is the velocity coefficient that controls the influence of the previous velocity, which may either be kept constant or let it vary in the next time steps: where iter is the current iteration number and itermax is the maximum number of iterations set for the algorithm run.According to this equation, kv decreases linearly to zero while k a increases to 2α as the number of iterations increases, which preserves the balance between the exploration and the speed of convergence [41].It is noteworthy that the parameters α and β in Equation ( 7) are tunable and defining these parameters result in definition of the acceleration and velocity coefficients (k a and kv).Value of 0.5 for both parameters α and β has been recommended in the reference paper of the CSS algorithm [41].
Substituting for k a and kv from Equation (7), Equations ( 5) and ( 6) can be rewritten as: In addition, to save the best results, a memory, known as the Charged Memory (CM), is recommended [41].If each CP moves out of the search space, its position is corrected using the harmony search-based handling approach, in which a new value is produced or selected from the CM, on a probabilistic basis.It is highly recommended to refer to the main reference [41] for better understanding the concepts and structure of the CSS algorithm, as some concepts in the present paper might not be described as detailed.
In the original CSS algorithm, when the calculations of the amount of forces are completed for all CPs, the new locations of the CPs are determined and also CM updating is fulfilled.In other words, the new location for each CP is determined after completion of an iteration and before commencement of the new iteration.Kaveh and Talatahari [45], ignoring this assumption, proposed the Enhanced CSS algorithm in which after evaluation of each CP, all updating processes are performed.Using this method of updating in the CSS algorithm, the new position of each agent can affect the moving process of the subsequent CPs while in the standard CSS, unless an iteration is completed, the new positions are not utilized.This enhanced algorithm, compared to the original CSS, while not requiring additional computational time, improves the performance of the algorithm by using the information obtained by CPs instantly.In a detailed investigation, considering the i th CP in the original CSS, although the solutions obtained by the CPs with a number less than i are created before the selected agent is used, however, these new designs cannot be employed to direct the i th CP in the current iteration.On the other hand, the original CSS archives the information obtained by the agents until a pre-determined time and this results in a break in the optimization process, while in the enhanced CSS algorithm the information of the new position of each agent is utilized in the subsequent search process, and this procedure improves the optimization abilities of the algorithm and also increases the convergence speed [45].

Water-Supply and Hydropower Reservoir Operation
In a water-supply reservoir operation, the objective is to obtain a set of releases from the reservoir (or a set of reservoir storage volumes) for the operation period with given inflow such that a predefined pattern of demands is met.In the other words, the objective is to set the released flow as close as possible to the demand and decrease the unnecessary overflows from the reservoirs, and hence, minimize the water deficit.Therefore, optimal operation of a water supply reservoir can be stated mathematically as [34,37]: 2 (10) subject to continuity equations at each time step: S min ≤ S(t) ≤ S max (12) where  (15) to the existing data.In cases where the evaporation loss is not considered in the formulations, Loss (t) is excluded from the Equation (11).
In a hydropower reservoir operation, the objective is to obtain a set of releases from the reservoir (or a set of reservoir storage volumes) such that the power generation from the reservoir is maximum, or as close as possible to the installed capacity of the hydro-electric plant.Hydropower operation of a single reservoir may be defined as [34,46]: subject to the continuity constraints defined by equations 11 to 15 defined for the simple operation problem.Here p (t) is power generated in megawatts (MW) in time step t, power is the installed capacity of hydro-electric plant (MW), and other parameters are defined as before.The power generated in time step t can be stated as follow: in which h (t) is the effective head of the hydroelectric plant as defined by Equation (18): H (t) is the elevation of water in reservoir at time step t which may be defined as a function of storage in the reservoir as: where g is the Earth's gravity acceleration, η is the efficiency of the hydroelectric plant, r (t) is release from reservoir (m 3 /s), PF is the plant factor, TWL (tail water level) is the downstream elevation of the hydroelectric plant (m), a, b, c and d are constants that can be obtained by fitting Equation (19) to the reservoir's data.

Model Performance Evaluation: Mathematical Test Functions
The CSS algorithm has previously been tested with some mathematical tests and structural optimization problems in the original reference [41].However, water resources management problems are more complex than structural problems and require a higher level of robustness in utilized approaches.In order to assess the performance of the proposed algorithm in complex and nonlinear problems in the field of water resources management, it is first applied to benchmark constrained and unconstrained mathematical optimization functions.For evaluation of robustness of the standard CSS algorithm and also assess the impact of utilizing the new method of updating on performance of the algorithm (called Enhanced CSS or ECSS in this section) [45], the functions are optimized with both standard CSS and ECSS.
Although value of 0.5 for both tunable coefficients in Equation ( 7) (i.e., α and β) have been recommended in the reference paper of the CSS algorithm [41], different values may result in better solutions.A sensitivity analysis for mathematical test functions shows that the value of 0.8 for these coefficients results in better solutions and increases the convergence rate of the algorithm.The best, worst, mean and standard deviation of the results are reported for 10 different runs of the algorithm for each of the mathematical test functions.
The first optimization problem is the Ackley's function [47], a continuous and multi-modal function defined by modulating an exponential function with a cosine wave of moderate amplitude.Ackley's function is defined as: where c 1 = 20, c 2 = 0.2, c 3 = 2π and n is taken equal to 2 here.The functions surface is a nearly flat outer region with moderate fluctuations converging to a hole in the middle.Multiple hills and valleys on the surface cause moderate complexity for optimization methods, as the search algorithms performing based on hill-climbing techniques are most likely to be trapped in local optima (Figure 1).An algorithm with a large scanning span that searches a wider neighborhood would be able to avoid the valleys and located better optima.Therefore, Ackley's function provides one of the reasonable test cases for the CSS algorithm.The results obtained by the CSS and ECSS algorithms as well as those obtained by Genetic Algorithm (GA) [48], Honey-Bees Mating Optimization algorithm (HBMO) [49] and Particle Swarm Optimization algorithm (PSO) [50] are presented in Table 1.The results for CSS and ECSS shown in Table 1 are best out of ten runs of the algorithms.The ECSS algorithm with 10 CPs, reaches the fitness value of 2 × 10 6 after 440 function evaluations (44 iterations with 10 CPs) and the best fitness value of 0 (global optimum) after 930 function evaluations in the best run of the algorithm.
Results show that all ten executions of the algorithms reach quite close to the global optimum value of the objective function, where the standard deviation of the objective function value over ten runs is 1.6 × 10 −15 for the CSS and 1.4 × 10 −15 for the ECSS.Table 1 demonstrates the impact of using the new method of updating in the structure of the algorithm on the convergence rate of the algorithm, as ECSS shows higher convergence speed compared to the standard CSS.The table denotes that the CSS algorithm obtains more accurate values in smaller number function evaluations in comparison with the other metaheuristic rivals, in terms of accuracy and convergence speed.The second numerical example is an unconstrained sine function defined as [48]: As seen in Figure 2a, the search space for this function is a highly non-linear and multi-modal surface.The ECSS algorithm with 30 CPs reached the best fitness value of 38.85029 after 1590 function evaluations (53 iterations with 30 CPs).The results obtained by the CSS and ECSS algorithms as well as those obtained by GA [48], HBMO [49] and PSO [50] are presented in Table 2.The results for CSS and ECSS shown in Table 2 are best out of ten runs of the algorithms.Results show that the standard deviation of objective function value is approximately zero, indicating that all 10 runs have converged to approximately one single solution.As seen from the Table 2, the CSS algorithm performs considerably faster comparing to the other rival approaches, where the CSS and ECSS algorithms locate the optima with 1950 and 1590 times function evaluations, respectively.
To assess the performance of the CSS algorithm in handling constrained problems, a twovariable, two-constraint constrained exponential function is considered [49], defined as (Figure 2b):  The second numerical example is an unconstrained sine function defined as [48]: As seen in Figure 2a, the search space for this function is a highly non-linear and multi-modal surface.The ECSS algorithm with 30 CPs reached the best fitness value of 38.85029 after 1590 function evaluations (53 iterations with 30 CPs).The results obtained by the CSS and ECSS algorithms as well as those obtained by GA [48], HBMO [49] and PSO [50] are presented in Table 2.The results for CSS and ECSS shown in Table 2 are best out of ten runs of the algorithms.Results show that the standard deviation of objective function value is approximately zero, indicating that all 10 runs have converged to approximately one single solution.As seen from the Table 2, the CSS algorithm performs considerably faster comparing to the other rival approaches, where the CSS and ECSS algorithms locate the optima with 1950 and 1590 times function evaluations, respectively.
To assess the performance of the CSS algorithm in handling constrained problems, a two-variable, two-constraint constrained exponential function is considered [49], defined as (Figure 2b): Subject to: The unconstrained objective function f (x1, x2) has a minimum solution at (3, 2) with a function value equal to zero.However, due to multiple constraints imposed to the function, this solution is not feasible and the constrained optimal solution is x = (2.2461,2.3815) with a function value equal to f = 13.61227.The feasible region is only approximately 0.7% of the total search space, which is a narrow crescent-shaped region, and the optimum solution is lying on the second constraint.
Employing the ECSS algorithm with 20 CPs, the best obtained fitness value is 13.59087 at x = (2.23809,2.24677)after 600 function evaluations (30 iterations with 20 CPs).Results obtained by the CSS and ECSS algorithms as well as those obtained by GA [48], HBMO [49] and PSO [50] are presented in Table 3.The results for CSS and ECSS shown in Table 3 are the best out of ten runs of the algorithms.All ten runs show a very small discrepancy with the global result as indicated by a very small value of the standard deviation.However, the ECSS reaches the global optima with less function evaluations.As seen in Table 3, the ECSS algorithm finds the optima with considerably less function evaluations comparing the other metaheuristic approaches.
Hydrology 2019, 6, x FOR PEER REVIEW 8 of 17 The unconstrained objective function f (x1, x2) has a minimum solution at (3, 2) with a function value equal to zero.However, due to multiple constraints imposed to the function, this solution is not feasible and the constrained optimal solution is x = (2.2461,2.3815) with a function value equal to f = 13.61227.The feasible region is only approximately 0.7% of the total search space, which is a narrow crescent-shaped region, and the optimum solution is lying on the second constraint.
Employing the ECSS algorithm with 20 CPs, the best obtained fitness value is 13.59087 at x = (2.23809,2.24677)after 600 function evaluations (30 iterations with 20 CPs).Results obtained by the CSS and ECSS algorithms as well as those obtained by GA [48], HBMO [49] and PSO [50] are presented in Table 3.The results for CSS and ECSS shown in Table 3 are the best out of ten runs of the algorithms.All ten runs show a very small discrepancy with the global result as indicated by a very small value of the standard deviation.However, the ECSS reaches the global optima with less function evaluations.As seen in Table 3, the ECSS algorithm finds the optima with considerably less function evaluations comparing the other metaheuristic approaches.The last test function used for investigation of performance of the CSS algorithm facing highly non-linear multi-variable problems is the well-known multimodal and continuous Fletcher-Powell function [47], which is a non-separable, non-linear, and irregular function, described as:  The last test function used for investigation of performance of the CSS algorithm facing highly non-linear multi-variable problems is the well-known multimodal and continuous Fletcher-Powell function [47], which is a non-separable, non-linear, and irregular function, described as: where x j is the decision variable, α j is a random coefficient within the range of -π and π, a ij and b ij are random coefficients within the range of −100 and 100, and n is the dimension of the function (Figure 3).The optimum point of the function is at x j = α j where the objection function value will be equal to zero.A 30-variable Fletcher-Powell function is chosen here to be optimized utilizing the CSS algorithm.The ECSS algorithm with 20 CPs, reaches the fitness value of 1246.37 after 2E5 function evaluations and the best fitness value of 440.29 after 2E6 function evaluations in the best run of the algorithm.The results obtained by ECSS algorithm as well as those obtained by the GA [48], HBMO [49], PSO [50] and Nonlinear Programing (NLP) using LINGO 8.0 software [49] are presented in Table 4.The ECSS algorithm obtains more optimum values in smaller number function evaluations in comparison with the other metaheuristic rivals.However, due to complexity of the Fletcher-Powell function, the algorithm does not converge to optimal values at every attempt, as shown by a large standard deviation value of 8805.65, for 10 runs of the algorithm.Multiple runs of the algorithm are needed to obtain the best optima achievable by the approach, as also seen in other algorithms [49].(29) ) ,..., 2 , 1 ( where xj is the decision variable, αj is a random coefficient within the range of -π and π, aij and bij are random coefficients within the range of -100 and 100, and n is the dimension of the function (Figure 3).The optimum point of the function is at xj = αj where the objection function value will be equal to zero.A 30-variable Fletcher-Powell function is chosen here to be optimized utilizing the CSS algorithm.The ECSS algorithm with 20 CPs, reaches the fitness value of 1246.37 after 2E5 function evaluations and the best fitness value of 440.29 after 2E6 function evaluations in the best run of the algorithm.The results obtained by ECSS algorithm as well as those obtained by the GA [48], HBMO [49], PSO [50] and Nonlinear Programing (NLP) using LINGO 8.0 software [49] are presented in Table 4.The ECSS algorithm obtains more optimum values in smaller number function evaluations in comparison with the other metaheuristic rivals.However, due to complexity of the Fletcher-Powell function, the algorithm does not converge to optimal values at every attempt, as shown by a fairly large standard deviation value of 8805.65, for 10 runs of the algorithm.Multiple runs of the algorithm are needed to obtain the best optima achievable by the approach, as also seen in other algorithms [49].Results demonstrate that the standard CSS algorithm without any improvements is a robust and fast convergent approach which outperformed its other metaheuristic rivals in optimization of complex and multimodal mathematical functions and seems to be capable of handling the highly nonlinear and non-convex problem of large-scale reservoir operation.Results also indicate that utilizing the new method of updating in the structure of algorithm, referred here as enhanced CSS or ECSS, improves its convergence speed.The remainder of this paper presents application of the new approach in water supply and hydropower reservoir operation problems.In all reservoir operation problems presented in the next sections, the enhanced CSS algorithm is utilized and will simply be called the CSS algorithm.

Reservoir Operation Case Study
To evaluate the performance of the CSS algorithm in solution of large-scale water resources management problems, water-supply and hydropower operation of "Dez" reservoir in southern Iran has been considered.Total storage capacity of "Dez" reservoir in pre-defined normal water level is 2510 MCM and the average inflow of the reservoir over 40 years period from 1970-2010 is 5900 MCM.The initial storage of the reservoir is taken equal to 1430 MCM.The maximum and minimum allowable storage volumes are set equal to 3340 and 830 MCM, respectively.The maximum and minimum monthly water release set equal to 1000 MCM and zero, respectively.The coefficients of the volume-elevation curve defined by Equation ( 19) are used as: a = 249.83364,b = 0.58720, c = −1.37 × 10 −5 and d = 1.526 × 10 −9 .The total installed capacity of hydroelectric power plant of the Dez reservoir is 650 MW, being operated with plant factor of 0.417 and 90% efficiency.The tail water level in downstream is assumed constant at 172 m above sea level.
These problems are solved here using the CSS algorithm for optimal monthly operation over 5, 20 and 40-year time spans, which are 60, 240 and 480 monthly periods, respectively.The water-supply operation and hydropower operation problems are solved using the CSS algorithm, separately, and results are shown for each operation type.The parameters of acceleration (ka) and velocity (kv) coefficients in Equation ( 7) are taken as 0.3, 0.3 and 0.5 for the 60, 240 and 480-month operational period problems, respectively, obtained via a tuning procedure.Tests show that for best results in reservoir operation optimization utilizing the CSS algorithm, the charged sphere radius (a) in the Equation ( 2) should be set close to zero (1E9).All the results presented here are obtained using 40, 100 and 1000 CPs for the 60, 240 and 480-month problems, respectively.The number of objective function evaluations is limited to 400,000 for all executions of the algorithm.
The water-supply and hydropower operation of "Dez" reservoir is first solved disregarding the effect of evaporation from the reservoir utilizing the CSS.The results of 5 executions of the algorithm are presented in Table 5. Disregarding evaporation, the NLP solver (LINGO 9.0) produces the objective function values of 20.6 and 45.4 for the hydropower operation over 240, and 480 periods [37].Disregarding evaporation, the value of 45.8 for 480 months hydropower operation of "Dez" reservoir, using the LINGO 9.0, was also reported by earlier studies [34].The CSS approach results in optimal solution of 21.4936 and 50.3394 for 240 and 480 period hydropower problems, respectively.The water-supply and hydropower operation of "Dez" reservoir is then solved considering the evaporation losses from the surface on the storage volume of the reservoir (Table 6).The inclusion of evaporation further increases the non-linearity of the model, in particular for the hydropower operation model for which the well-known LINGO NLP solver has failed to find a feasible solution for 240 and 480 operational periods, as previously reported [46].The LINGO NLP solver (LINGO 9.0) has yielded optimum solutions of 0.732, 4.77, and 10.50 for the water-supply problem for 60, 240, and 480 months [46].The LINGO NLP solver was able to locate a near-optimal solution of 7.37 for the hydropower operation over the shortest operation period of 60 months, but failed to yield any feasible solution for the longer operation periods of 240 and 480 months, which may be due to the non-convexity of the hydropower operation, which is even higher for the 240 and 480-month problem [46].Failure in finding a feasible solution for 480-period hydropower operation of "Dez" reservoir considering the evaporation from the reservoir, using the NLP method, was reported by other studies as well [34].
Gradient-based nonlinear programming (NLP) methods solve problems that do not involve high level of nonlinearity in objective function and constraints.However, in models with large number of decision variables and/or highly nonlinearity, these approaches tend to fail in locating feasible solutions, or converge to local optima [51].In long-period reservoir operations, or in models with evaporation losses consideration, the non-linearity and non-convexity of the model rises and gradient-based NLP solvers may not be a suitable choice since they may either produce local suboptimal solutions and/or may even fail to locate any feasible solutions.
It can be seen from Table 6 that in all short and long operation period of 60, 240 and 480 months, the CSS algorithm can yield the near-optimal solutions in both water-supply and hydropower operation models.In longer operation periods of 240 and 480 months for the water-supply operation model, the CSS algorithm results in objective function values of 2.7302 and 9.4263 respectively.Comparing the results of the CSS algorithm with those obtained by LINGO 9.0 NLP solver demonstrates that the CSS approach outperforms the NLP method with the solutions of 4.77 and 10.5 for 240 and 480-month operation periods, respectively.In 240 and 480-month hydropower operation, the CSS results in near optimal solutions with objective function values of 22.9670 and 54.3472, respectively.The NLP method fails in producing any feasible solution for these problems, as noted earlier.These results indicate robustness of the CSS algorithm to solve both convex and non-convex, small-scale and large-scale, reservoir-operation problems.Table 7 presents the results obtained by genetic algorithm (GA), Particle Swarm Optimization (PSO), Ant Colony Optimization (ACO) and Big Bang-Big Crunch (BB-BC) algorithm for the current Dez reservoir operation problem [37,46,52].Comparing the results of Tables 6 and 7 shows that in all operation problems, these algorithm yield solutions inferior to those obtained by the CSS approach, with the difference being more significant for the longer operation periods.For the 480-month operations, the solutions yielded by the GA and PSO algorithms are significantly far from the optimal solution.
All 10 runs of the CSS algorithm for simple and hydropower operations have resulted in feasible solutions; while some of the other metaheuristic approaches failed in locating feasible solutions for all executions.These results can be compared with those obtained by the conventional Ant Colony Optimization Algorithm (ACO) [52].The results show that ACO was capable of producing 10 feasible solutions for the simplest case of water-supply operation over 60 monthly periods and 8 feasible solutions for the hydropower operation over 60 monthly periods.In longer operation periods, i.e., 240 and 480 monthly periods, the number of runs with a feasible solution decrease.For 240 monthly periods, only 8 and 7 feasible solutions were created for water-supply and hydropower operations, respectively while for 480 monthly periods, ACO was only capable of producing one feasible solution for both water-supply and hydropower operation.Table 7 shows that while GA was unable to find a feasible solution for the longest operation period, the PSO algorithm could produce feasible solution only for the shortest operation period [37,52].also show the ability to locate feasible solutions in all 10 runs of the algorithms, with inferior results compared to the CSS algorithm.It should be noted that the number of objective function evaluations was limited to 400000 for the CBB-BC algorithm as well [37].Figure 4a shows the best solution obtained by the CSS algorithm (monthly releases) for 60-month water-supply operation, considering evaporation losses, versus the monthly water demand defined by the problem.Figure 4b illustrates the storage volume of the reservoir at each time step, calculated based on the releases ruled by the CSS algorithm and given river inflow.The Figure 4 shows the storage volume at each time step is confined between maximum (S_max) and minimum (S_min) allowable storage, defined by the problem.

Conclusions
Water resources management problems, including optimal reservoir operation, are complex and nonlinear optimization problems with large number of decision variables, which require search for more robust methods with high convergence speed for their solution.In this study, a robust metaheuristic optimization algorithm called the Charged System Search algorithm was introduced to the field of water resources management, and was applied to a large-scale real-world reservoir operation problem for the first time.The efficiency of the CSS algorithm was pre-evaluated using four well defined and highly nonlinear benchmark mathematical functions.Results show that the CSS approach outperforms its metaheuristic rivals considered in this study, in terms of accuracy and convergence speed.The CSS was applied to the optimization of water-supply and hydropower operation of the "Dez" reservoir in Iran for operation periods of 60, 240, and 480 months, and the results compared with those derived by Genetic, PSO, ACO, and CBB-BC algorithms as well as those obtained by the NLP approach.The results suggest that the CSS algorithm solves large-scale reservoir-operation problems with higher accuracy and convergence speed compared to the other available heuristic-search methods considered in this study.The results also indicated superiority of the CSS algorithm to the LINGO 9.0 NLP solver in highly nonlinear problem of extracting operation policies for the long-term evaporation-included hydropower operation model of "Dez" reservoir, in which the NLP solver fails to yield a feasible solution.As the non-linearity of the problem was increased by adding evaporation and/or expanding the operation period, the CSS approach outperforms the NLP solver.The CSS algorithm proves to be a robust, accurate, and fast convergent approach in handling complex water resources problems, as the effectiveness of this algorithm was shown in the present study.However, performance of this algorithm in other water management problems, such as multi-reservoir operation problems and conjunctive surface/ground water resources management, is yet to be investigated.Additionally, the reservoir operation problem presented in this study was formulated based on known river inflow values to the reservoir.However, uncertainties in river flow may affect the water release policies.Inclusion of such uncertainties in reservoir operation formulations may further improve the results to achieve more reliable operation rules.

Figure 2 .
Figure 2. Surface of the unconstrained sine function (a) and the constrained exponential function (b).

Figure 2 .
Figure 2. Surface of the unconstrained sine function (a) and the constrained exponential function (b).

Figure 4 .
Figure 4. Best solution yielded by the CSS algorithm for 60-month water-supply operation considering evaporation losses: water demand versus release from the reservoir (a), and river inflow versus storage volume at each time step (b).As seen from the figure, the storage is confined between maximum (S_max) and minimum (S_min) allowable storage.

Figure 4 . 17 Figure 5 .
Figure 4. Best solution yielded by the CSS algorithm for 60-month water-supply operation considering evaporation losses: water demand versus release from the reservoir (a), and river inflow versus storage volume at each time step (b).As seen from the figure, the storage is confined between maximum (S_max) and minimum (S_min) allowable storage.Hydrology 2019, 6, x FOR PEER REVIEW 14 of 17

Figure 5
Figure5presents variation of the objective function value versus the number of function evaluations for the best solution obtained by the CSS algorithm for "Dez" reservoir water-supply and hydropower operation, over each operational period.

Figure 5 .
Figure 5. Convergence curve of the optimum solution obtained by the CSS algorithm for "Dez" reservoir water-supply and hydropower operation, considering evaporation losses, over 60 (a), 240 (b) and 480 (c) months periods.
is the number of time steps, D (t) is water demand in time step t in million cubic meters (MCM), R (t) is release from the reservoir in time step t (MCM), Dmax is maximum demand (MCM), S (t) is storage at the start of time step t (MCM), I (t) is inflow in time step t (MCM), Smin and Smax are minimum and maximum storage of reservoir (MCM), respectively and Rmin and Rmax are minimum and maximum allowed release from reservoir (MCM), respectively.Loss (t) is net amount of gain and loss of the reservoir resulting from precipitation and evaporation in time step t.Ev(t) is the evaporation height during the time step t, and x 0 , x 1 , x 2 and x 3 are constants that can be obtained by fitting Equation

Table 1 .
Results of the standard Charged System Search (CSS) and Enhanced Charged System Search (ECSS) algorithm, in comparison with the other rival metaheuristics algorithms for Ackley's function.

Table 2 .
Results of the standard CSS and ECSS algorithms in comparison with the other heuristics for the unconstrained sine function.

Table 3 .
Results of the CSS and ECSS algorithms in comparison with the other heuristics for the constrained exponential function.

Table 2 .
Results of the standard CSS and ECSS algorithms in comparison with the other heuristics for the unconstrained sine function.

Table 3 .
Results of the CSS and ECSS algorithms in comparison with the other heuristics for the constrained exponential function.

Table 4 .
Results of the ECSS algorithm in comparison with the other rival metaheuristics algorithms for 30-variable Fletcher-Powell function.

Table 5 .
Results of "Dez" reservoir operation using CSS algorithm over 5-Disregarding evaporation loss.

Table 6 .
Results of "Dez" reservoir operation using CSS algorithm over 10 runs-Considering evaporation loss.