Parameter estimation of Muskingum model using grey wolf optimizer algorithm

Highlights • The flood routing is carried out by two non-linear Muskingum model.• The main purpose of this work is to make a comprehensive study between models optimized by AGWO, GWO and other meta-heuristic algorithms.• In order to compare the results of the GWO algorithm to those of more recent algorithms, the flood routing was performed by using the Augmented Grey Wolf Optimizer algorithm as well.• The GWO method generated good results on application to four case studies.


Background
Flood routing possesses a number of applications, such as the evaluation of peak flood discharge, the design of hydraulic structures, flood control and management, the drop in property damage as well as in both mental and physical injuries, and the design of flood warning systems. The flood routing methods are divided into two major classes of hydraulic and hydrological methods; within the latter category, the Muskingum method is the most common approach used for flood routing in both linear and non-linear forms. The number of hydrological parameters depends on the type of the Muskingum model employed. Niazkar and Afzali [27] and Kang et al. [16] have proposed different types of the Muskingum model. In addition, various researchers have optimized the hydrological parameters of the Muskingum model through several independent techniques, while the effort s to find more precise methods are ongoing. In the meantime, the meta-heuristic algorithms have been extensively used in the optimal determination of the Muskingum model parameters to date, by means of which it has always been attempted to carry out the flood routing more accurately.

Literature review
In recent decades, the artificial intelligence-based optimization methods have been widely utilized in the parameter estimation of the Muskingum model, whose outputs have been more accurate than the results of the Lagrange multiplier (LMM), segmented least square (S-LSM) methods. Mohan [24] and Kim et al. [18] implemented the genetic algorithm (GA) and the harmony search (HS) algorithm, respectively; in the former, a new vector is originated from merely two vectors, whereas every single vector in the latter forms a new one. Using the particle swarm optimization (PSO) algorithm, Chu and Chang [6] showed that the PSO outputs are more suitable than those of the GA approach, while the HS algorithm offers more precise results, in comparison with the PSO algorithm. The immune clonal selection algorithm (ICSA) features a higher convergence speed than traditional algorithms, which partly compensates for the defects associated with the GA. In another study, Luo and Xie [19] compared the performance of this algorithm to that of the HS and GA algorithms, and found that the ICSA has a decent speed. Geem [11] and Karahan et al. [17] utilized the parameter setting free-harmony search (PSF-HS) and harmony search-Broyden-Fletcher-Goldfarb-Shanno (HS-BFGS) algorithms, which, owing to their high convergence speed, are used in complex systems. In order to estimate the parameters of their Muskingum model, Niazkar and Afzali [26] adopted the honey-bee mating optimization (HBMO) algorithm, and compared the outputs to the results of 17 other algorithms. Noteworthy in this algorithm is the fast convergence to the optimum value within a broad range of values. This algorithm was subsequently employed by Niazkar and Afzali [27] in conjunction with the generalized reduced gradient (GRG) method for a six-parameter Muskingum model. Hamedi [14] investigated the invasive weed optimization (IWO) algorithm. Moghaddam et al. [23] estimated the Muskingum model parameters by using the PSO algorithm. In order to determine the four hydrological parameters of their Muskingum model, Farahani et al. [9] used the Shark algorithm. Norouzi and Bazargan [ 28 , 29 ] employed the PSO algorithm to estimate the parameters of Muskingum model.
The GWO algorithm is straightforward programming-wise, which also has a proper accuracy in engineering problems. Investigating the current literature disclosed that the GWO algorithm has not yet been utilized to estimate the parameters of nonlinear Muskingum models. In various publications, this algorithm has either been adopted as single or hybrid with other metaheuristic algorithms so as to improve the precision of the metaheuristic algorithms. Benefiting from the GWO algorithm, Mustafa et al. [25] optimized the parameters of the least-squares support-vector machine (LSSVM), and estimated the water level via the GWO-LSSVM algorithm. In another research work, Marroufpoor et al. [21] used the ANFIS-GWO algorithm for the soil moisture prediction, which was then compared to the artificial neural network (ANN), support vector regression (SVR) and adaptive neuro-fuzzy inference system (ANFIS) algorithms. Hamour et al. [15] adopted the Augmented Grey Wolf Optimizer algorithm (AGWO) to minimize total power losses, and compared their results to those of SPSO (Selective Particle Swarm Optimization) and FEB (Fuzzy Adaptation of Evolutionary Programming). They found that AGWO is an effective and competent method for radial distribution system reconfiguration. Yue et al. [39] developed the GWO algorithm in combination with the FWA (Fireworks Algorithm), named FWGWO, utilized for 16 benchmark functions, whose results were then compared to those of nine other algorithms. Tikhamarine [33] used the GWO algorithm, along with an artificial neural network, for the monthly prediction of the stream flow.

Innovation and objectives
As mentioned in the previous section, various methods have been proposed over the recent decade to improve the Muskingum model's accuracy. For instance, the optimization algorithm can be employed to estimate the Muskingum model parameters. The use of various algorithms directly affects the objective function.The purpose of this paper is to utilize the AGWO and GWO algorithm in order to calculate the optimal hydrological parameters of two three-and four-parameter Muskingum models, along with comparing the results obtained to those of other algorithms applied in flood routing. The aforementioned algorithms have thus been employed in four case studies, in each of which the outflow hydrograph was developed after having determined the Muskingum model hydrological parameters through the AGWO and GWO algorithms, whose outputs were then compared to those of other algorithms. It has been indicated that this two algorithms determines the hydrological parameters of the three-and four-parameter Muskingum models more precisely, in comparison with all the other algorithms used. It should be noted that although it has outperformed the GWO algorithm in other publications, the results of the AGWO algorithm in all the four case studies are similar to those of the GWO algorithm, with no superiority over each other in terms of the accuracy of the results obtained. This is despite the fact that the GWO algorithm has reached its convergence point in less iteration number and population size, on account of which its outputs were employed in this paper.

Non-linear Muskingum model
The Muskingum model was originally proposed by the United States Army Corps of Engineers, in 1938, for the flood control in the area along the Muskingum river [1] . Equations (1) and (2) are the nonlinear Muskingum models respectively with three and four parameters (named NL3 and NL4, respectively): where S t = channel storage at time t (L 3 ); I t = inflow at time t (L 3 /T); O t = outflow at time t (L 3 /T); k = time-storage constant for river basin; x = a weighting parameter usually in the range of 0 to 0.3; while a and m = exponent parameters to consider the nonlinearity effect. Equation (1) can be obtained in case a = 1 in Eq. (2) . Assuming a = 1, once again, leads to the outflow magnitude in Eq. (1) . Using Eq. (2) : The continuity equation is as follows: Substituting Eq. (3) into Eq. (4) results in: The objective function, the minimizer of the sum of squared errors among the routed and observed ouflows ( SSQ ), can be computed as follows: where O t : observed outflow and ˆ O t : routed outflow. The flood routing steps are recommended in the following order: 1) A value is assigned to k, x, m and a ; 2) The initial value of S t is calculated from Eq. (1) ; 3) Changes in the storage over time (time derivative) is evaluated by using Eq. (5) ; 4) The storage capacity can then result from the following equation: 5) The outflow is then computed from Eq. (3) -in many research works, the terms ( I t + I t−1 2 ) or I t−1 has been used instead of I t ; 6) The SSQ value is calculated; 7) The magnitudes of x, k, m and a are updated based on the optimization algorithm; 8) Stages 2-7 should be reiterated until reaching the termination condition [ 3 , 11 ] The evaluation indices of SAD, EQ p , ET p , MARE and VarexQ have been used as per the following equations in order to assess the performance of the algorithm in this study, as well as to compare it to the previously mentioned methods.
In these equations, O P and ˆ O P represent the peak of observed and routed outflows; T P and ˆ T P denote the peak time at observed and calculated peak outflows, respectively; and Ō t is mean of observed outflows. The less the EQ p is, the more precise will be the model results, while the lower ET p lead to a more accurate estimation of the peak discharge occurred. Equation (22) can be applied to measure the shape proximity and size of hydrographs [ 2 , 23 ].

Grey wolf optimization (GWO) algorithm
The grey wolf algorithm, originally developed by Mirjalili et al. [22] , is a subset of swarm intelligence algorithms, which has a hierarchical structure. There are four ranks in each grey wolf pack, the first of which belongs to the alpha male wolves that assert dominance on other groups; the second category comprises the beta male wolves, which enhance the decisions of the first group through giving feedbacks to the alpha male wolves; the next group includes the omega wolves. If a wolf is not an alpha, beta or omega, it is called delta. Delta wolves are supposed to submit to alphas and betas, yet they dominate the omega. The grey wolf algorithm has been established on the basis of pack hunting, which is of interest to grey wolves. The hunting process encompasses three main phases: a) tracking, chasing, and approaching the prey; b) pursuing, surrounding, and exhausting the prey until it stops moving; and c) attacking the prey.
Initially, the population decision variables are randomly developed from the search agents, which are basically the grey wolves, by considering the upper and lower bounds, following which the population of the wolves are sorted from lower to higher fitness based on the value of the fitness function ( SSQ ). The first member in the sorted population is called the alpha wolf, the second one is the beta wolf, and the third one is the delta wolf. The distance of the alpha, beta and delta wolves from the prey is calculated via Eqs. (13 -16 ).
D δ , D β and D a are the distances of the alpha, beta and delta wolves from the prey; t is the number of iterations; C and X symbolize the vectors of coefficients and spatial location, respectively; and r 1i are random vectors within the range of [0,1]. The spatial location of the alpha, beta and delta wolves can be updated based on Eqs (17 -20 ).
Where A is the vector of coefficients, and r 2i represents the random vectors within the range of [0,1]. The components of the vector a decrease linearly from two to zero.
Different spatial locations close to the best option could be achieved by adjusting the vectors A and C . In order to simulate the hunting process mathematically, it is initially assumed that the alpha, beta and delta wolves have a better knowledge during the hunt. Therefore, the first three optimal solutions are saved, in accordance with which the other search agent wolves (the omegas) are obliged to update their spatial location. The new population members can be evaluated by averaging the new spatial locations as per the following equation: The best member of the new population is introduced as the solution. The process of developing a new population, along with choosing its best member as the solution, continues until meeting the termination condition [22]  The augmented grey wolf algorithm, developed by Qais et al. [32] , has augmented the exploration and exploitation process of grey wolf optimizer algorithm. In this algorithm, only alpha and beta wolves are used for hunting and the components of the vector a decrease nonlinearly from two to one according to Eq. (22) .

Results and discussions
In this research, there are four case studies on which the aforementioned method was tested: a) Wilson River; b) River Wye; c) Viessman and Lewis [36] ; and d) Karun River. The flood routing results obtained from the NL3 and NL4 models, as well as those optimized by the GWO algorithm, are as follows: Case study 1 -The Wilson data has been selected as the first example to compare the performance of the nonlinear Muskingum model in this research to those adopted in various studies [37] . With the time step of six hours, the hydrograph considered here has a single peak with an initial inflow rate of 22 cms (Cubic meter per second), which reaches the maximum value of 111 cms in a period of 30 hours, and then diminishes to 18 in 96 hours. Alternatively, the outflows were calculated after having estimated the hydrological parameters in the NL3 and NL4 models via the GWO algorithm. In Table 1 , the results of the NL3 and NL4 models are compared to those of twenty NL3 models and five NL4 models. In the former category, the optimal values of k, x and m as well as the objective function value ( SSQ ) have been calculated as 0.51745, 0.28689, 1.8681 and 36.7679, in that order. In the latter category, the optimal values of k, x, m and a are 0.1391, 0.2956, 4.0830 and 0.4326, respectively; and the SSQ value is equal to 7.6674. It should be noted that the time step of the NL3 and NL4 models was considered to be 6 hours and one hour, respectively. As evidenced in Table 1 Table 2 presents the values of the routed outflows for both models, through which the inflow and outflow hydrographs, together with the routed outflow hydrograph, have been plotted in Fig. 1 , in which the close resemblance between the hydrograph routed by means of the GWO algorithm and the observed outflow hydrograph is evident. Using the GWO algorithm, the values of MARE, VarexQ, EQ T , EQ p were optimized for the NL3 and NL4 models, Table 10 . Among the various methods applied in the optimal estimation of the hydrological parameters of the NL3 and NL4 models, the GWO algorithm has well reduced the value of MARE and ET p ; and the VarexQ value is at its maximum, which indicates the good fitness of the proposed model  to the observed outflows. Considering the EQ P value, it can also be declared that the GWO algorithm, similar to some other ones, is well capable of estimating the peak discharge.
Case study 2 -This is a non-smooth flood hydrograph belonging to the River Wye in the United Kingdom, which has been studied in many research works. Table 3 ( and GRG algorithms, as well as by 12%, 9%, 9% and 8% for the NL4 model, compared to the GA-GRG, SFLA-NMS, LINGO and PSO algorithms, respectively. Increasing the number of hydrological parameters from three to four has reduced the SSQ value by 11%. For the NL3 and NL4 models, the values of the routed outflows and their corresponding hydrograph can be found in Table 4 and Fig. 2 , respectively. Moreover, the performance evaluation indices of the algorithms are shown in Table 10 , according to which the SAD, MARE and EQ p parameters have less values than those in other algorithms, whereas the VarexQ value has increased. The EQ p was obtained 0.07, being the lowest value among other methods and indicating that the peak discharge was estimated with lower error compared to other methods. Case study 3 -A multi-peak flood hydrograph has been used in this example [36] . Table 5 compares the GWO algorithm to the AGWO, IICSA, GA-GRG, PSO, SFLA-NMS and LINGO algorithms. In the NL3 model, the optimal values of the hydrological parameters estimated by GWO, which reduced the SSQ value by 29% than the IICSA algorithm, are k = 0.4562 , x = 0.4250 , m = 1.2540 and SSQ = 51742 , respectively. These were obtained as k = 0. 4535 , x = 0.4266 , m = 1.1987, a = 1.0470 and SSQ = 51527 for the NL4 model. Comparing the current GWO outputs to those of the previous research discloses that this algorithm has diminished the SSQ value by 33%, 31%, 30% and 30% with respect to the GA-GRG, PSO, SFLA-NMS and LINGO algorithms. Table 6 presents the routed outflow values and the optimal hydrological parameters for the NL3 and NL4 models. In addition, the hydrographs of the routed and observed outflows and inflows are depicted in Fig. 3 . Table 10 shows the evaluation indices calculated via the GWO and five other algorithms. Noting the values of SAD, EQ p , ET p , MARE and VarexQ in this example, it can be claimed that GWO estimates the outflows and their peak more efficiently  than the other aforementioned algorithms. The VarexQ is also higher in this algorithm, which indicates the better agreement between the routed and observed outflows. Thus, the GWO algorithm in this example offers more suitable results in comparison with those in the previous studies.
Case study 4 -After having been proven in achieving the desired outcomes in the previous examples, the GWO algorithm was adopted for a real river in Iran. Karun is the longest river in the southwest of Iran, to which the Dez river -another important river in the ecosystem -connects as well. The target discharges were recorded by two gauging stations called Godar and Gotvand, with the minimum and maximum inflows being 380 and 1300 cms, respectively. In this single peak hydrograph, the time step is equal to two hours [23] , while this was considered to be one hour in the routing process. The parameter estimation results via the GWO, AGWO and GA, ABC, SA, SFLA and PSO algorithms are tabulated in Table 7 , in which the optimal value of the hydrological parameters of the NL3 model have been obtained as k = 14982.7, x = 0.1457 and m = 0.1361, with their corresponding SSQ as 59294. In this model, the GWO algorithm outperformed the GA, ABC, SA and SFLA algorithms, in terms of substantially reducing the objective function value to the extent that it has decreased by 68%, 67%, 56% and 55%, compared to the mentioned algorithms, respectively. The significant decrease in the objective function is a strong point of the GWO algorithm when compared to a highly capable algorithm such as GA. In the NL4 model, the SSQ value of the GWO algorithm declined by 18% than that of the PSO algorithm; and the associated optimal hydrological parameters were evaluated as k = 41735.3, x = 0.0762, m = 0.0210 and a = 3.393, with the corresponding SSQ value of 56698. In the NL3 and NL4 models, in order to assess and further compare the GWO method to other evolutionary algorithms, the performance criteria indices presented in Table 10 were calculated by using the optimal values of the hydrological parameters attained through the GWO algorithm. As evidenced in Table 10 , the GWO method has a superior performance; As shown in Table 10 , the values of VarexQ and SAD were obtained 98.39 and 1077, respectively. The values of these indexes show that the GWO algorithm had better performance compared to the other methods listed in the table. The inflow and outflow hydrographs of Karun river and the routed outflow hydrograph, as well as the routed outflow values to provide a step-by-step and more accurate comparison, are presented, respectively, in Fig. 4 and Table 8 , in which there is an insignificant difference between the routed and observed outflows. It should be added that the model has been fitted well to the observational data.
As evidenced in Tables 1-8 , similar to the GWO algorithm, the AGWO algorithm has estimated the hydrological parameters in all the four case studies with a good accuracy as well, the results of which are also superior to those of the other algorithms presented in Tables 1-8 .  Table 9 compares these two algorithms in terms of their population size and maximum iteration number. Therefore, the final outputs in this research have been reported based on the GWO algorithm.

Conclusion
Given the significance of flood forecasting and its role in designing the flood-resistant hydraulic structures, as well as the measures required for flood control and management, the flood routing is of special interest to the researchers in this field, to the extent that it has always been attempted to achieve this imperative quickly and accurately. In this research, the flood routing was conducted by using two nonlinear Muskingum models with three and four constant parameters. In order to estimate the parameters of both models, the grey wolf optimization (GWO) algorithm was utilized in optimization problems because of its efficiency, which has not yet been applied in flood routing processes; and it was thus used for three case studies and a real example in Iran (Karun River), in each of which its performance has been compared to that of other evolutionary algorithms, in terms of the SSQ, SAD, EQ p , ET p , MARE and VarexQ values. It was revealed, in all the four case studies, that the GWO algorithm is well capable of optimally estimating the Muskingum model parameters, whose application in flood routing problems is strongly recommended. Eventually, a comparison between the results of the AGWO and GWO algorithms has been conducted in four case studies. The aforesaid algorithms in the four case studies have reduced the SSQ value up to 55% than the minimum SSQ value, which can be implemented to optimize the hydrological parameters of other forms of the Muskingum equation, whether single or combined with other algorithms such as GA and PSO.