Abstract
In this study, the capability of the recently introduced moth swarm algorithm (MSA) was compared with two robust metaheuristic algorithms: the harmony search (HS) algorithm and the imperialist competitive algorithm (ICA). First, the performance of these algorithms was assessed by seven benchmark functions having 2–30 dimensions. Next, they were compared for optimization of the complex problem of fourreservoir and 10reservoir systems operation. Furthermore, the results of these algorithms were compared with nine other metaheuristic algorithms. Sensitivity analysis was performed to determine the appropriate values of the algorithms’ parameters. The statistical indices coefficient of determination (R^{2}), root mean square error (RMSE), mean absolute error (MAE), mean square error (MSE), normalized MSE (NMSE), mean absolute percentage error (MAPE), and Willmott’s index of agreement (d) were used to compare the algorithms’ performance. The results showed that MSA was the superior algorithm for solving all benchmark functions in terms of obtaining the optimal value and saving CPU usage. ICA and HS were ranked next. When the dimensions of the problem were increased, the performance of ICA and HS dropped but MSA has still performed extremely well. In addition, the minimum CPU usage and the best solutions for the optimal operation of the fourreservoir system were obtained by MSA, with values of 269.7 seconds and 308.83, which are very close to the global optimum solution. Corresponding values for ICA were 486.73 seconds and 306.47 and for HS were 638.61 seconds and 264.61, which ranked them next. Similar results were observed for the 10reservoir system; the CPU time and optimal value obtained by MSA were 722.5 seconds and 1,195.58 while for ICA they were 1,421.62 seconds and 1,136.22 and for HS they were 1,963.41 seconds and 1,060.76. The R^{2} and RMSE values achieved by MSA were 0.951 and 0.528 for the fourreservoir system and 0.985 and 0.521 for the 10reservoir system, which demonstrated the outstanding performance of this algorithm in the optimal operation of multireservoir systems. In a general comparison, it was concluded that among the 12 algorithms investigated, MSA was the best, and it is recommended as a robust and promising tool in the optimal operation of multireservoir systems.
HIGHLIGHTS
This is the first application of robust MSA in the optimization of water resources.
The capability of 11 metaheuristic algorithms in optimization of multireservoir systems was compared.
MSA was the superior model in terms of CPU time and obtaining the optimal solution.
INTRODUCTION
The optimal operation of reservoirs is one the most important issues in water resources management, especially for multireservoir systems. Many control variables determine the operation strategies for scheduling a sequence of releases to meet a large number of demands for different users. Optimal reservoir operation defines the optimum policies to retain or release water from a reservoir at different periods of the year according to the inflows and demands.
To develop a comprehensive and optimal operation policy for a reservoir, the sources of uncertainties should be considered. The uncertainties may be due to different factors such as inaccurate predictive models, data scarcity, measurement and observation errors, the measurement site being unrepresentative, or problems in aggregating or disaggregating data. For a multireservoir system, the uncertainties may be caused by the stochastic nature of the system inputs (rainfall, reservoir inflow, evaporation, leakage, release, etc.), the nonlinearity of functions, the multiple constraints, the large number of decision variables, and other spatial and temporal variations of system components. Consideration of these uncertainties can help in the design of efficient operational policies and to develop robust predictive models.
Over the past few years, several evolutionary algorithms have been developed and applied to solving reservoir optimization problems (Table 1).
Evolutionary algorithm .  References .  Problem definition .  Results . 

Genetic algorithm (GA)  Sharif & Wardlaw (2000) 


Chang et al. (2005) 

 
GarousiNejad et al. (2016) 

 
Genetic programming (GP)  FallahMehdipour et al. (2013a, 2013b) 


FallahMehdipour et al. (2012) 

 
Ashofteh et al. (2015) 

 
FallahMehdipour et al. (2013a, 2013b) 

 
Ant colony optimization (ACO)  Jalali et al. (2006) 


Kumar & Reddy (2006) 

 
Moeini & Afshar (2013) 

 
Mohammed et al. (2018) 

 
Particle swarm optimization (PSO)  Afshar (2013) 


Ghimire & Reddy (2013) 

 
Afshar (2012) 

 
AlAqeeli & Agha (2020) 

 
Harmony search (HS)  BashiriAtrabi et al. (2015) 


Kougias & Theodossiou (2013) 

 
Mirbeyk et al. (2020) 

 
Water cycle algorithm (WCA)  Bozorg Haddad et al. (2015) 


Qaderi et al. (2018) 

 
Honeybee mating optimization (HBMO)  Bozorg Haddad et al. (2011) 


Soghrati & Moeini (2020) 

 
Gravity search algorithm (GSA)  BozorgHaddad et al. (2016) 


Imperialist competitive algorithm (ICA)  Afshar et al. (2015) 


HosseiniMoghari et al. (2015) 

 
Rahimi et al. (2020) 

 
Spider monkey algorithm (SMA)  Ehteram et al. (2018) 


Bat algorithm (BA)  Ahmadianfar et al. (2016) 


Evolutionary algorithm .  References .  Problem definition .  Results . 

Genetic algorithm (GA)  Sharif & Wardlaw (2000) 


Chang et al. (2005) 

 
GarousiNejad et al. (2016) 

 
Genetic programming (GP)  FallahMehdipour et al. (2013a, 2013b) 


FallahMehdipour et al. (2012) 

 
Ashofteh et al. (2015) 

 
FallahMehdipour et al. (2013a, 2013b) 

 
Ant colony optimization (ACO)  Jalali et al. (2006) 


Kumar & Reddy (2006) 

 
Moeini & Afshar (2013) 

 
Mohammed et al. (2018) 

 
Particle swarm optimization (PSO)  Afshar (2013) 


Ghimire & Reddy (2013) 

 
Afshar (2012) 

 
AlAqeeli & Agha (2020) 

 
Harmony search (HS)  BashiriAtrabi et al. (2015) 


Kougias & Theodossiou (2013) 

 
Mirbeyk et al. (2020) 

 
Water cycle algorithm (WCA)  Bozorg Haddad et al. (2015) 


Qaderi et al. (2018) 

 
Honeybee mating optimization (HBMO)  Bozorg Haddad et al. (2011) 


Soghrati & Moeini (2020) 

 
Gravity search algorithm (GSA)  BozorgHaddad et al. (2016) 


Imperialist competitive algorithm (ICA)  Afshar et al. (2015) 


HosseiniMoghari et al. (2015) 

 
Rahimi et al. (2020) 

 
Spider monkey algorithm (SMA)  Ehteram et al. (2018) 


Bat algorithm (BA)  Ahmadianfar et al. (2016) 


As can be seen, several evolutionary algorithms have been developed in recent years. Because of the many advantages of these algorithms, such as rapid rate of convergence, the capability of a large number of simulations to arrive at an optimum solution, avoiding local optima, the ability for multiobjective optimization, easy handling of nonlinearity and nonconvexity of the problem domain, lack of restriction by the number of dimensions and computational requirements, they have become an attractive alternative to the classical methods for the optimization of engineering problems.
The main question that this paper seeks to answer is that, among the several evolutionary algorithms that have been recently introduced by different researchers in different engineering fields, which one is the most appropriate for the optimal operation of multireservoir systems? Accordingly, this paper investigates the capability of the recently introduced moth swarm algorithm (MSA) for the optimal operation of four and 10reservoir systems in comparison to the harmony search (HS) algorithm and the imperialist competitive algorithm (ICA). To the authors' knowledge, this is the first application of MSA for the optimal operation of multireservoir systems. The extraordinary performance of these three algorithms, as outlined in Table 2, has made them the most successful algorithms among dozens of algorithms in solving and optimizing various complex engineering problems (Yang et al. 2017; Yoo et al. 2018; Gerist & Maheri 2019; Kim et al. 2019; Afsari et al. 2020; Khamari et al. 2020; Sayed et al. 2020). The results of these three algorithms will be also compared with nine other algorithms developed by previous researchers, including GA, honeybee mating optimization (HBMO), gravity search algorithm (GSA), improved bat algorithm (IBA), FA, modified FA (MFA), krill, hybrid GAkrill, WCA.
Algorithm .  Advantages .  Limitations . 

MSA 


ICA 


HS 


Algorithm .  Advantages .  Limitations . 

MSA 


ICA 


HS 


MATERIALS AND METHODS
As previously stated, this study investigates the capability of three robust algorithms, MSA, HS and ICA, for the optimal operation of multireservoir systems. The overall methodology of this study is shown in Figure 1. More details of methodology are presented below.
Moth swarm algorithm (MSA)
MSA is a robust metaheuristic algorithm which originates from moth behaviour in nature. Moths hide from predators during the day, and at night they use celestial navigation to orient themselves in the dark and exploit food sources. In celestial navigation, the flying direction lies at a constant angle to the moonlight as a remote light source. In this algorithm, the best solution is indicated by the position of the light source, and the quality of this solution is determined by the intensity of the light source. Each moth swarm includes three groups: pathfinders, prospectors and onlookers. Pathfinders can explore the best point over the optimization space with the firstinlastout principle to lead the movement of the main swarm. Prospectors wander into a random path close to the light sources, which have been marked by the pathfinders. Onlookers move directly toward the best global solution (moonlight), which has been obtained by prospectors.
The optimization procedure in MSA is such that, at any iteration, each moth is incorporated into the optimization problem to find the highest intensity of its corresponding light source. The best fitness in the swarm is considered as the pathfinders’ positions, and guidance for the next updated iteration. Hence, the second and third best groups are the prospectors and onlookers, respectively. The optimization of the algorithm is performed through three phases of initialization, reconnaissance and transverse orientation. In initialization, at the beginning of the flight, the initial positions of moths are randomly generated. Next, the moths are grouped into the three types (pathfinders, prospectors, onlookers) based on their calculated fitness. To prevent a probable premature convergence in this phase as well as to improve the diversity of solutions, pathfinders update their positions by crossover operations and Levy mutation. During this procedure, the type of each moth changes dynamically. For instance, if a prospector finds a solution that is brighter than the existing one, it is promoted to become a pathfinder moth. During the optimization procedures, by decreasing the number of prospectors, the number of onlookers increases, which may lead to a global solution. In this phase, MSA forces the onlookers to search more effectively (by celestial navigation) and find the optimal solution. The governing equations and mathematical explanation of MSA was provided by Mohamed et al. (2017). Less than 3 years after development of this algorithm, its extraordinary performance has been reported by several researchers (Zhou et al. 2018; Shilaja & Arunprasath 2019; Duong et al. 2020; Kotb & ElFergany 2020).
Harmony search (HS) algorithm
HS is an optimization algorithm inspired by the musical improvization process in music bands. It was first introduced by Geem et al. (2001), and in recent years, has been widely employed as an optimization tool in various engineering problems. Harmony in music is similar to the optimization solution vector and the musicians’ improvisations whereby the musicians’ search for harmony is analogous to local and global search schemes in the optimization procedure. HS has many advantages (e.g. fewer mathematical requirements, no need for initial value setting of the decision variables, having stochastic random searches, and strong mechanism to generate new solution vectors) which made it a robust and prominent algorithm in the optimization of complex engineering problems. It has been successfully employed in the optimization of truss structural optimizations (Cheng et al. 2016; Saka et al. 2016), engineering management optimizations (Shen et al. 2017), economic load dispatch problem (AlBetar et al. 2016); clustering applications (Abualigah et al. 2020; Talaei et al. 2020), electrical engineering (Shiva & Kumar 2020), medical sciences (Koti et al. 2020), civil and geotechnical engineering (Bekdas et al. 2020), optimal design of water distribution networks (Geem 2006; Geem & Cho 2011), optimal operation of multireservoir systems (Geem 2007), simulation of irrigation systems (Cisty 2008), optimization of groundwater management (Ayvaz 2009), identification of unknown groundwater pollution sources (Ayvaz 2010) and optimal operation of reservoir system with emphasis on flood forecasting (BashiriAtrabi et al. 2015).
Imperialist competitive algorithm (ICA)
ICA was proposed by AtashpazGargari & Lucas (2007). It is based on a global heuristic search that applies imperialism and imperialistic competition process as a source of inspiration. ICA starts with an initial population. Some of the best individuals of this population, called countries, are picked up as the imperialist states and all the rest become the colonies of these imperialists. Due to the imperialists' powers, which are inversely proportional to their cost, the colonies of the initial population are divided among them. Having distributed the colonies between the imperialists and established the initial empires, these colonies commence proceeding toward their relevant imperialist country. After the successful application of ICA in different engineering optimizations, it has been approved by several researchers as one of the most powerful optimization algorithms (AbdElazim & Ali 2016; Al Dossary & Nasrabadi 2016; Fathy & Rezk 2017; Mikaeil et al. 2018; Emami & Parsa 2020; Lei et al. 2020).
All the algorithms were coded in MATLAB 2016a using a PC with i7 CPU 1.8 GHz/16GB RAM/2TB HDD.
Verification of the algorithms
In order to evaluate the efficiency and validation of the developed models for the optimal operation of multireservoir systems, a set of standard benchmark functions was selected, as presented in Table 3.
Benchmark function .  Dimensions .  Global optimum .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
GoldstenPrice  2  3  3  5.51  3  9.23  3  5.13 
McCormick  2  −1.9132  −1.9132  4.97  −1.9132  8.36  −1.9132  6.92 
Rosenbrock  2  0  0  5.53  1.8 × 10^{−6}  12.68  2.7 × 10^{−15}  7.14 
Rosenbrock  10  0  0  5.92  1.43  16.46  0.055  10.24 
Rosenbrock  30  0  2.66 × 10^{−6}  6.06  63.03  21.24  2.11  11.34 
Shekel  4  −10.5364  −10.5321  5.53  −6.4619  11.62  −10.532  8.21 
Sixhump camel  2  −1.0316  −1.0316  4.94  −1.0316  9.91  −1.0316  5.13 
Benchmark function .  Dimensions .  Global optimum .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
GoldstenPrice  2  3  3  5.51  3  9.23  3  5.13 
McCormick  2  −1.9132  −1.9132  4.97  −1.9132  8.36  −1.9132  6.92 
Rosenbrock  2  0  0  5.53  1.8 × 10^{−6}  12.68  2.7 × 10^{−15}  7.14 
Rosenbrock  10  0  0  5.92  1.43  16.46  0.055  10.24 
Rosenbrock  30  0  2.66 × 10^{−6}  6.06  63.03  21.24  2.11  11.34 
Shekel  4  −10.5364  −10.5321  5.53  −6.4619  11.62  −10.532  8.21 
Sixhump camel  2  −1.0316  −1.0316  4.94  −1.0316  9.91  −1.0316  5.13 
The performance of MSA in solving these functions was compared with HS and ICA. The population size and the number of evaluations of benchmark functions in all the algorithms was identical and proportional to the dimensions of each function. As seen in Table 3, for benchmark functions with lower dimensions, the performance of all the algorithms was approximately similar, but for large dimension problems (Rosenbrock function with dimension of 10 and 30), MSA was the only algorithm which was capable of solving the problem. The performance of the other algorithms was so weak, the results of them dramatically diverged from the optimal value. Figure 2 shows the convergence rate of the algorithms in reaching the optimal value for standard benchmark functions. It indicates that for all benchmark functions, MSA reached its optimal value by the smallest number of iterations. The difference was more obvious in high dimension problems.
After successful verification of utilized algorithms by benchmark functions, their performance was investigated in the optimization of two multireservoir systems.
Sensitivity analysis of model parameters
Before using the algorithms to optimize the multireservoir system operation, a sensitivity analysis was performed to determine the best values for the algorithm parameters. The tuning parameters include number of iterations, number of variables, number of search agents and number of pathfinders (for MSA), number of iterations, number of variables, harmony memory size, harmony memory considering rate and pitchadjusting rate (for HS), and number of iterations, number of variables, number of countries, revolution rate and assimilation coefficient (for ICA). The sensitivity of each model to all corresponding parameters was investigated and the best values for model parameters were identified by sensitivity analysis.
Test cases: optimization of multireservoir systems
The fourreservoir problem was first formulated and solved by Larson (1968) as a linear problem with a known global optimum. This case was later used by a number of researchers for the assessment of different optimization algorithms (Wardlaw & Sharif 1999; HosseiniMoghari et al. 2015; Ehteram et al. 2017). The fourreservoir system, as illustrated in Figure 3(a), consists of both series and parallel connections. The supplies from the system are used for hydropower generation and irrigation. Hydropower generation is possible from each reservoir. The outflow from reservoir four may be diverted for irrigation. Hydropower and irrigation benefits are quantified by linear functions of discharge. There are inflows to the first and second reservoirs only, and these are 2 and 3 units, respectively, in all time periods. The initial storage in all reservoirs is 5 units.
The 10reservoir problem is more complicated, not only in terms of size, but also because of the many timedependent constraints on storage. The layout of the 10reservoir problem is shown in Figure 3(b). The system comprises reservoirs in series and in parallel, and a reservoir may receive supplies from one or more upstream reservoirs. Inflows are defined for each of the upstream reservoirs, and initial storage and target storage at the end of the operating period are specified for each reservoir. In addition, there are minimum operating storage amounts in each reservoir that must be satisfied, as well as constraints on minimum and maximum reservoir releases. Operation of the system is optimized over 12 operating periods to maximize hydropower production. More details of these problems can be found in Murray & Yakowitz (1979).
Here S_{k}(t) = storage at time t in reservoir k (k = 1, …., K); I_{k}(t) = inflows in time period t to reservoir k; = minimum storage in reservoir k; = maximum storage in reservoir k; = minimum release from reservoir k; = maximum release from reservoir k.
Evaluation criteria
In the above equations, is the releases in time period t from the investigated algorithms, is the optimum release in time period t, is the mean of the optimum release, is the mean of the release, and n is the number of total time periods. The variation domain of Willmott's index of agreement ranges from to 1, so that 1 indicates perfect agreement between the optimum release and releases from the investigated algorithms.
A low value of RMSE and a high value of R^{2} indicate acceptable accuracy of the algorithm and correlation between the data, and also imply its superiority over the other algorithms. Each of the MAE, MSE, NMSE, and MAPE parameters shows the difference between the optimum release and releases from the investigated algorithms; the lower the values of these parameters, the more efficienct the algorithm. MSE highlights the difference between the data, and its normalized form (NMSE) can be compared with other algorithms.
RESULTS AND DISCUSSION
Results of sensitivity analysis
As previously described, to ensure the reliability and validity of the models’ outputs, MSA, HS and ICA were put through 10 runs with different iterations. To achieve the appropriate number of iterations, a sensitivity analysis on the number of iterations was performed. Accordingly, each algorithm was run with 500, 1,000, 2,000 and 5,000 iterations. Table 4 shows the values of objective function for different numbers of iterations in the fourreservoir system.
Iteration .  Algorithm .  

MSA .  HS .  ICA .  
500  307.5855  244.4548  281.3062 
1,000  308.25  258.842  296.363 
2,000  308.094  259.1158  298.869 
5,000  308.731  261.4383  299.71 
Iteration .  Algorithm .  

MSA .  HS .  ICA .  
500  307.5855  244.4548  281.3062 
1,000  308.25  258.842  296.363 
2,000  308.094  259.1158  298.869 
5,000  308.731  261.4383  299.71 
Only MSA could reach the values very close to the optimal value after 1,000 iterations. Although the two other algorithms converged after 1,000 iterations, they have a large difference from the optimal value. For iterations over 1,000, while the execution time of the algorithms dramatically increased, the variations of the objective function were negligible. Therefore, the number of iterations in each algorithm was considered to be 1,000. The results of sensitivity analysis for the fourreservoir system are shown in Figures S1–S4 in the Supplementary Information. In addition, Figure 4 shows the variations of objective function by increasing the number of iterations. Similar results were obtained for the 10reservoir system. Table 5 shows the final values of algorithm parameters for multireservoir systems, all of which were derived by sensitivity analysis.
MSA .  Parameter .  Number of iterations .  Number of variables .  Number of search agents .  Number of pathfinders .  – . 

Value .  1,000 .  48–120 .  100 .  20 .  – .  
.  Parameter .  Number of iterations .  Number of variables .  Harmony memory size .  Harmony memory considering rate .  Pitchadjusting rate . 
HS .  Value .  1,000 .  48–120 .  100 .  0.95 .  0.3 . 
.  Parameter .  Number of iterations .  Number of variables .  Number of countries .  Revolution rate .  Assimilation coefficient . 
ICA .  Value .  1,000 .  48–120 .  100 .  0.4 .  0.5 . 
MSA .  Parameter .  Number of iterations .  Number of variables .  Number of search agents .  Number of pathfinders .  – . 

Value .  1,000 .  48–120 .  100 .  20 .  – .  
.  Parameter .  Number of iterations .  Number of variables .  Harmony memory size .  Harmony memory considering rate .  Pitchadjusting rate . 
HS .  Value .  1,000 .  48–120 .  100 .  0.95 .  0.3 . 
.  Parameter .  Number of iterations .  Number of variables .  Number of countries .  Revolution rate .  Assimilation coefficient . 
ICA .  Value .  1,000 .  48–120 .  100 .  0.4 .  0.5 . 
Convergence rate and models performance
Figure 5 shows the convergence rate of applied algorithms to reach the optimum value in the operation of four and 10reservoir systems. It indicates the rapid convergence of MSA compared to HS and ICA. The results in Figure 5 are consistent with those in Figure 2, in which MSA is the best and HS is the worst algorithm for solving the largescale problems. It can be seen from these from the figures that MSA reaches the optimum value in fewer iterations, but HS has premature convergence. Premature convergence means the solution is stuck within the local optimum.
Tables 6 and 7 show the values of objective function and the average CPU run time of each algorithm for 10 runs of the algorithms for the fourreservoir and 10reservoir systems, respectively. It was found that MSA was able to produce the best solutions for both the four and 10reservoir systems. The results of 10 different runs of MSA generated an optimum solution close to 100% of the global optimum for both multireservoir systems. As can be seen in Table 6, the value of the objective function achieved by MSA was 16.03% and 3.86% better than that of HS and ICA, respectively. The calculation time was 638 seconds for HS, 486 seconds for ICA, and 269 seconds for MSA, indicating the outstanding speed of MSA execution. Thus, MSA ranks best in terms of the objective function production and saving calculation time. Also, according to Table 7, the value of the objective function for the 10reservoir system achieved by MSA was 12.5% and 6.6% better than that of HS and ICA, respectively. The calculation time of the 10reservoir system optimization by MSA, HS and ICA were 722 seconds, 1,963 seconds, and 1,421 seconds, respectively, indicating the superior performance of MSA.
Number of runs .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
1  308.83  376.12  246.36  678.67  304.14  524.12 
2  308.52  337.13  256.74  751.64  292.67  498.26 
3  308.64  447.01  261.01  716.23  295.12  543.56 
4  307.52  447.90  261.41  649.32  293.99  551.72 
5  307.70  383.85  252.74  675.92  292.06  486.73 
6  308.47  308.41  264.61  706.64  306.47  601.24 
7  308.43  269.71  258.51  668.65  296.53  524.93 
8  307.41  313.91  262.84  645.34  294.58  512.65 
9  308.83  318.64  260.35  638.61  293.8  503.15 
10  308.14  318.54  263.85  706.54  294.27  561.23 
Best  308.83  264.61  306.47  
Worst  307.41  246.36  292.06  
Average  308.25  258.842  296.363  
SD  0.53  5.631914614  4.900340238  
Coefficient of variation  0.0017  0.021758117  0.016534926  
Best CPU time (s)  269.71  638.61  486.73 
Number of runs .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
1  308.83  376.12  246.36  678.67  304.14  524.12 
2  308.52  337.13  256.74  751.64  292.67  498.26 
3  308.64  447.01  261.01  716.23  295.12  543.56 
4  307.52  447.90  261.41  649.32  293.99  551.72 
5  307.70  383.85  252.74  675.92  292.06  486.73 
6  308.47  308.41  264.61  706.64  306.47  601.24 
7  308.43  269.71  258.51  668.65  296.53  524.93 
8  307.41  313.91  262.84  645.34  294.58  512.65 
9  308.83  318.64  260.35  638.61  293.8  503.15 
10  308.14  318.54  263.85  706.54  294.27  561.23 
Best  308.83  264.61  306.47  
Worst  307.41  246.36  292.06  
Average  308.25  258.842  296.363  
SD  0.53  5.631914614  4.900340238  
Coefficient of variation  0.0017  0.021758117  0.016534926  
Best CPU time (s)  269.71  638.61  486.73 
Number of runs .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
1  1,195.13  827.51  1,024.86  2,132.64  1,125.12  1,532.74 
2  1,195.56  812.12  1,034.24  2,036.45  1,091.23  1,462.15 
3  1,194.12  723.34  1,060.76  1,963.41  1,103.41  1,421.62 
4  1,192.10  722.55  1,035.67  2,009.35  1,109.37  1,574.62 
5  1,193.36  914.79  1,060.41  2,103.62  1,110.65  1,596.38 
6  1,192.77  983.29  1,048.24  1,992.72  1,136.22  1,496.21 
7  1,194.93  851.16  1,053.09  2,012.04  1,113.84  1,466.56 
8  1,195.58  826.97  1,040.36  2,031.26  1,107.39  1,514.81 
9  1,192.78  751.16  1,032.84  1,984.32  1,132.72  1,441.26 
10  1,195.56  739.63  1,054.78  2,101.36  1,121.63  1,526.31 
Best  1,195.58  1,060.76  1,136.22  
Worst  1,192.10  1,024.86  1,091.23  
Average  1,194.19  1,044.525  1,115.158  
Standard deviation  1.34  12.61843823  13.82279341  
Coefficient of variation  0.0011  0.012080552  0.012395368  
Best CPU time (s)  722.55  1,963.41  1,421.62 
Number of runs .  MSA .  HS .  ICA .  

Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  Optimal value .  CPU time (s) .  
1  1,195.13  827.51  1,024.86  2,132.64  1,125.12  1,532.74 
2  1,195.56  812.12  1,034.24  2,036.45  1,091.23  1,462.15 
3  1,194.12  723.34  1,060.76  1,963.41  1,103.41  1,421.62 
4  1,192.10  722.55  1,035.67  2,009.35  1,109.37  1,574.62 
5  1,193.36  914.79  1,060.41  2,103.62  1,110.65  1,596.38 
6  1,192.77  983.29  1,048.24  1,992.72  1,136.22  1,496.21 
7  1,194.93  851.16  1,053.09  2,012.04  1,113.84  1,466.56 
8  1,195.58  826.97  1,040.36  2,031.26  1,107.39  1,514.81 
9  1,192.78  751.16  1,032.84  1,984.32  1,132.72  1,441.26 
10  1,195.56  739.63  1,054.78  2,101.36  1,121.63  1,526.31 
Best  1,195.58  1,060.76  1,136.22  
Worst  1,192.10  1,024.86  1,091.23  
Average  1,194.19  1,044.525  1,115.158  
Standard deviation  1.34  12.61843823  13.82279341  
Coefficient of variation  0.0011  0.012080552  0.012395368  
Best CPU time (s)  722.55  1,963.41  1,421.62 
As previously described, for a better comparison of the three utilized algorithms in the optimal operation of multireservoir systems, seven statistical indices were employed. According to Table 8, the maximum values of the accuracy parameters R^{2} and d were obtained by MSA, with values of 0.9511 and 0.9874 for the fourreservoir system and 0.9852 and 0.9963 for the 10reservoir system. The accuracy indices dropped significantly in HS and ICA for both multireservoir systems. In the same way, while the minimum values of error parameters were achieved by MSA in both the fourreservoir (RMSE = 0.5275, MAE = 0.1951, MSE = 0.3506, NMSE = 0.04993 and MAPE = 490.1845) and the 10reservoir (RMSE = 0.5211, MAE = 0.1357, MSE = 0.4254, NMSE = 0.0150 and MAPE = 282.0782) systems, the corresponding values for HS and ICA were not satisfactory.
Case study .  Algorithm .  Performance metrics .  

R^{2} .  RMSE .  MAE .  d .  MSE .  NMSE .  MAPE .  
Fourreservoir system  MSA  0.9511  0.5275  0.1951  0.9874  0.3506  0.04993  490.1845 
HS  0.5475  1.7826  1.3346  0.8442  3.1776  0.4525  9,489.9  
ICA  0.8108  1.1838  0.5288  0.9493  1.4015  0.1996  6,259.4729  
10reservoir system  MSA  0.9852  0.5211  0.1357  0.9963  0.4254  0.0150  282.0782 
HS  0.4151  4.1406  2.7930  0.7933  17.1447  0.6043  17,878  
ICA  0.6824  3.0784  1.9390  0.9087  9.4765  0.334  12,073 
Case study .  Algorithm .  Performance metrics .  

R^{2} .  RMSE .  MAE .  d .  MSE .  NMSE .  MAPE .  
Fourreservoir system  MSA  0.9511  0.5275  0.1951  0.9874  0.3506  0.04993  490.1845 
HS  0.5475  1.7826  1.3346  0.8442  3.1776  0.4525  9,489.9  
ICA  0.8108  1.1838  0.5288  0.9493  1.4015  0.1996  6,259.4729  
10reservoir system  MSA  0.9852  0.5211  0.1357  0.9963  0.4254  0.0150  282.0782 
HS  0.4151  4.1406  2.7930  0.7933  17.1447  0.6043  17,878  
ICA  0.6824  3.0784  1.9390  0.9087  9.4765  0.334  12,073 
Figure 6(a) and 6(b) show the water release patterns for the operation of four and 10reservoir systems using MSA, HS and ICA. As can be seen, the operating policies obtained using MSA gave the maximum benefits with a more appropriate release pattern, so that the system does not face shortages. Regarding the higher capability of MSA in calculating nearoptimal global solutions, it can be employed by water policymakers as a guide (rule curve) to schedule water releases from multireservoir systems in a way that gives the most benefits. It can be seen from Figure 6 that HS and ICA failed to produce reasonable results for multireservoir systems.
Comparison with previous research
Table 9 compares the results of this study with those of previous studies on the optimization of four and 10reservoir system operation. It can be seen that MSA has the highest performance among all evolutionary algorithms. For the fourreservoir system, MSA gave the optimal value of 308.83, the closest to the global optimum. After that, GSA (BozorgHaddad et al. 2016), with the optimal value of 308.30, demonstrated better performance than the other algorithms. Next were MFA (GarousiNejad et al. 2016), hybrid GAkrill (Ehteram et al. 2017), IBA (Ahmadianfar et al. 2016), HBMO (BozorgHaddad et al. 2011), krill (Ehteram et al. 2017), ICA (present study), WCA (Eskandar et al. 2012; BozorgHaddad et al. 2015; Qader et al. 2018), FA (GarousiNejad et al. 2016), and HS (present study) with the optimum values of 308.21, 308.17, 308.05, 307.50, 307.26, 306.47, 306.39, 305.51 and 264.61, respectively. This confirms the superiority of MSA to the 11 other metaheuristic algorithms in the optimal operation of the fourreservoir system.
Study .  Algorithm .  Average value of the objective functions .  

Fourreservoir system .  10reservoir system .  
Wardlaw & Sharif (1999)  GA  –  1,190.25 
BozorgHaddad et al. (2011, 2016)  HBMO  307.50  1,148.05 
GSA  308.30  –  
Ahmadianfar et al. (2016)  IBA  308.05  1,192.89 
GarousiNejad et al. (2016)  FA  305.51  1,097.41 
MFA  308.21  1,183.59  
Ehteram et al. (2017)  Krill  307.26  1,189.66 
Hybrid GAkrill  308.17  1,193.91  
Qaderi et al. (2018)  WCA  306.39  1,172.42 
Present study  HS  264.61  1,060.76 
ICA  306.47  1,136.22  
MSA  308.83  1,195.58 
Study .  Algorithm .  Average value of the objective functions .  

Fourreservoir system .  10reservoir system .  
Wardlaw & Sharif (1999)  GA  –  1,190.25 
BozorgHaddad et al. (2011, 2016)  HBMO  307.50  1,148.05 
GSA  308.30  –  
Ahmadianfar et al. (2016)  IBA  308.05  1,192.89 
GarousiNejad et al. (2016)  FA  305.51  1,097.41 
MFA  308.21  1,183.59  
Ehteram et al. (2017)  Krill  307.26  1,189.66 
Hybrid GAkrill  308.17  1,193.91  
Qaderi et al. (2018)  WCA  306.39  1,172.42 
Present study  HS  264.61  1,060.76 
ICA  306.47  1,136.22  
MSA  308.83  1,195.58 
For the case of the 10reservoir system, it was also observed that MSA was the most robust model. It achieved the optimal value of 1,195.58, close to 100% of the global optimum. After MSA, the GAkrill hybrid (Ehteram et al. 2017) solved the problem with the highest percentage of the objective function relative to the global optimal solution. The order of next algorithms in terms of achieving the best solutions IBA, GA, krill, MFA, WCA, HBMO, ICA, FA, and HS, respectively. Comparison of the results indicates the extraordinary capability of MSA in the optimization of a complex 10reservoir system. It could efficiently reach the absolute global optimum in both multireservoir systems.
CONCLUSION
In this study, the capability of the recently introduced MSA for solving benchmark functions as well as the optimal operation of multireservoir systems was compared with two robust algorithms: ICA and HS. The results from MSA, HS and ICA were compared with those of nine other metaheuristic algorithms: GA, HBMO, GSA, IBA, FA, MFA, krill, hybrid GAkrill and WCA, which were developed by previous researchers. The results indicated the remarkable performance of MSA compared to all other algorithms in the optimal operation of fourreservoir and 10reservoir systems. This algorithm successfully reached the absolute global optimum in both multireservoir systems. It also had the shortest CPU time for obtaining the optimal value. The findings of this research and the application of the same methodology for the optimal operation of multireservoir systems will enable decision makers to make informed choices on water development, conservation, allocation, and use in the context of growing demands for all uses and increased scarcity. In addition to the low cost, easy implementation and simple procedure of the methodology presented, it has many capabilities that make it attractive to be used by water policymakers, water resources planners, and reservoir managers.
ACKNOWLEDGEMENTS
We are grateful to the Research Council of Shahid Chamran University of Ahvaz for financial support (GN: SCU.WH98.26878). Also, we acknowledge the Khuzestan Water and Power Authority.
DATA AVAILABILITY STATEMENT
All relevant data are included in the paper or its Supplementary Information.