Multi-Objective Parameter Estimation of Improved Muskingum Model by Wolf Pack Algorithm and Its Application in Upper Hanjiang River , China

In order to overcome the problems in the parameter estimation of the Muskingum model, this paper introduces a new swarm intelligence optimization algorithm—Wolf Pack Algorithm (WPA). A new multi-objective function is designed by considering the weighted sum of absolute difference (SAD) and determination coefficient of the flood process. The WPA, its solving steps of calibration, and the model parameters are designed emphatically based on the basic principle of the algorithm. The performance of this algorithm is compared to the Trial Algorithm (TA) and Particle Swarm Optimization (PSO). Results of the application of these approaches with actual data from the downstream of Ankang River in Hanjiang River indicate that the WPA has a higher precision than other techniques and, thus, the WPA is an efficient alternative technique to estimate the parameters of the Muskingum model. The research results provide a new method for the parameter estimation of the Muskingum model, which is of great practical significance to improving the accuracy of river channel flood routing.


Introduction
In recent years, frequent floods have caused serious losses to people's lives and property [1].Therefore, it is urgent to identify the flood routing rules to reduce flood losses.As a classic solution of flood routing, the Muskingum model continues to be a popular method for flood routing [2].The Muskingum model was first developed by McCarthy [3].In practical application, the key problem for applying the Muskingum model is parameter estimation [4,5], which is a highly nonlinear optimization problem.The studies on the Muskingum model mainly include the parameter estimation method and model improvement.
Traditional parameter estimation methods mainly include the Trial Algorithm (TA), the moment method, the least squares method, and the differential algorithm [5,6].However, these methods are limited by the optimal estimation of the channel storage curve, which has led to a significant error between the calculated results and the observed data.In recent years, intelligence algorithms have been widely used in solving the nonlinear problem accurately and efficiently [7].Many researchers have applied various techniques to estimate the parameters of the Muskingum model in recent years.For example, Luo and Xie [4] applied the immune clonal selection algorithm to estimate the parameters of the Muskingum model, getting a higher precision than the other techniques.Mohon [8] proposed the genetic algorithm to estimate the parameters of two nonlinear Muskingum Water 2018, 10, 1415 2 of 14 routing models.Chen and Yang [9] applied the Gray-encoded accelerating genetic algorithm in the parameters optimization of the Muskingum model.Chu and Chang [10] used the Particle Swarm Optimization (PSO) algorithm to estimate the Muskingum model parameters, which improves the accuracy of the Muskingum model for flood routing.Barati [11] applied the parameter-setting-free technique, interfaced with a harmony search algorithm to the parameter estimation of the Muskingum model, and found good model parameter values.Zhang et al. [12] applied the Shuffle Complex Evolution Algorithm (SCE-UA) to optimize and estimate the discharge proportion coefficient x of the Muskingum equation and the number of streams subsections partitioned for mainstreams and tributaries.Niazkar and Afzali [13] proposed a hybrid method, which combined the Modified Honey Bee Mating Optimization and Generalized Reduced Gradient algorithms, reducing the sum of the squared (SSQ) value for the double-peak case study.Ouyang et al. [14] applied the invasive weed optimization (HIWO) algorithm to the parameter estimation of nonlinear Muskingum models.Although various techniques were applied to estimate the parameters of the nonlinear Muskingum flood routing model, the application of the model still suffers from a lack of an efficient method for parameter estimation.An efficient method for parameter estimation in the Muskingum model is still required to improve computational precision.
The improved Muskingum model mainly includes parameter setting and optimization criteria [15,16], and studies on model improvement have had a significant breakthrough in recent years.Zhang et al. and Vatankhah [17,18], for example, used a nonlinear Muskingum flood routing model with variable exponent parameters, producing the most accurate fit for outflow data.Moghaddam et al. [19] proposed a new Muskingum model with four parameters, the sum of the squared (SSQ) or absolute (SAD) deviations between the observed and estimated outflows considered as objective functions.Although the new Muskingum model becomes more complex, it improves the fit to observed flow, especially in multiple-peak hydrographs.Luo et al. [20] proposed a multi-objective estimation routine of the Muskingum model, involving single-peak, multi-peak, and non-smooth hydrographs, proving that the multi-objective estimation procedure is consistent and effective in estimating the parameters of the Muskingum model.Easa [21] pointed out models that adopt the outflow criterion result in a poor fit to the observed storages, presenting a new approach that incorporates both criteria in the estimation process and aids trade-off analysis.
This paper proposes a parameter estimation of the Muskingum model based on a new intelligent algorithm-the Wolf Pack Algorithm (WPA).In this paper, several floods with different magnitudes from river channels (Ankang hydropower station-Ankang city and Ankang hydropower station-Shuhe) of the Hanjiang River were first selected to estimate their respective parameters.Secondly, an improved multi-objective of the Muskingum model as well as that of the WPA and its solving steps are proposed.Thirdly, the results from the different algorithms in terms of the weighted sum of absolute difference and the coefficient of determination between the observed and routed outflows are compared to verify the performance of the WPA.Finally, conclusions are drawn based on the results.

Research Area
The Ankang to Shuhe section of Hanjiang River in China was chosen as the study area (Figure 1).The length of the river section is 109.36km, accounting for 7% of the total reach of the Han River, with a drop of 61 m and an average gradient ratio of 0.06% of the river bed.Due to the steep slope of the mountain basin, the poor permeability of the rock formation, and poor adjustment of the channel, the flood process is very fast, the peak shape is thin, the inter-annual variability of the flood is larger, and the flood variation coefficient is larger.The parameters of the river and flood data are shown in Tables 1 and 2.   The Muskingum model is a traditional method for solving river flood routing, which is mainly solved by the continuity equation and the dynamic equation [22].In order to obtain the flood routing equation, the water balance equation and storage equation are solved by:  The Muskingum model is a traditional method for solving river flood routing, which is mainly solved by the continuity equation and the dynamic equation [22].In order to obtain the flood routing equation, the water balance equation and storage equation are solved by: where Q 1 and Q 2 are the outflow of the downstream section at the beginning and the end of the period, respectively, I 1 and I 2 are the inflow of upstream section at the beginning and the end of the period, K is channel storage coefficient, and x is specific gravity coefficient of flow.C 0 , C 1 , and C 2 are the flow routing coefficients.

Design of Objective Function
The optional parameters of the Muskingum model can be x and K, or C 0 , C 1 , and C 2 .If K and x are chosen as the optimization variables, it is necessary to ensure that the water storage capacity of the river is a single linear relationship with the reservoir flow during the optimization process, which is difficult to determine, given the range of variables.Meanwhile, there is no need to go deep into the physical parameters of the model, because of the intelligent algorithm based on the black box model [23].Therefore, C 0 and C 1 are chosen as the optimization variables.The results of river flood routing are mainly reflected in the flood process and the flood peak of the simulation.
The results of river flood flow evolution are mainly reflected in the degree of fitting of flood process and flood peak to the actual flood.Some studies indicated that the success of a calibration process is highly dependent on the objective function chosen as a calibration criterion [19].The most commonly used objective function for the calibration procedure is the SSQ errors between observed and computed outflow [4,5,20,24], but some research has indicated that the SSQ is not necessarily correct [16,25].Considering the above reasons, the objective function established in this paper can be described as follows: Objective 1: Minimum weighted sum of absolute difference (SAD): where Q 0 (i) is the actual outflow of the downstream section at a time i, I c (i − 1) and I c (i) are the simulated flow for the upstream section at time (i − 1) and i, respectively, Q 0 (i − 1) is the simulated flow for the upstream section at time (i − 1), and n is the length of the time series used for calibration.C 0 , C 1 , and C 2 are the flow routing coefficients.
The SAD will give the minimum difference between observed and computed outflow [9,19,26].Objective 1 is the SAD multiplying the corresponding weight taken from the observed flow at the corresponding time, which will increase large flow influence on the parameter estimation, especially on the flood peak.Thus, the weight can increase the simulation accuracy of the flood peak.
Objective 2: Maximum coefficient of determination D c : where D c is the coefficient of determination, Q c (i) is the computed value at time i, Q 0 (i) is the observed value at time i, Q 0 is the mean value of the observed flood process, and n is the number of time series used for calibration [22].The deterministic coefficient is an indicator to measure the consistency between the flood forecast and the observed process.In this paper, the coefficient of determination as Objective 2 can ensure that the simulated flow process is the closest to the observed values [27,28].
In order to reduce the multi-objective calculation and adaptation algorithm, the comprehensive multi-objective function chosen in this paper is as follows: where f is the comprehensive objective, f 1 and f 2 are Objective 1 and Objective 2, respectively, and φ 1 and φ 2 are the weight coefficient of different objectives, respectively.

Wolf Pack Algorithm (WPA)
The Wolf Pack Algorithm is a novel swarm intelligence algorithm with strong local search ability and global convergence.It consists of the leading wolf, fierce wolves, and explore wolves, including three kinds of intelligent behaviors: Scouting, summoning, and beleaguering.Simultaneously, with a productive rule for the leading wolf, which is that the winner can dominate all, it is a renewable mechanism, namely survival of the fittest, for a pack of wolves.It has been applied in the field of mathematics, physics, and hydropower station optimization, achieving good calculation results [29][30][31].In this paper, the WPA is used for the parameter estimation of the Muskingum model.Combining the principle of the WPA with the objective function, we designed a WPA for solving Muskingum model parameters.The procedure of our algorithm for parameter optimal estimation of the Muskingum model is shown as follows: Step 1: (Initial wolves) We must initially determine the size of the algorithm population n and the number of iterations gen, the probe wolf scale factor a, the population regeneration factor b, the step factor S, the distance determination factor W, and the maximum number of scouting T max .The Model Parameter (C 0 , C 1 ) is regarded as the position of the artificial wolf (x 1 , x 2 ) in a two-dimensional decision space.The position of the wolf pack is initialized in the range of C 0 and C 1 , as shown in Equation ( 7): where x i,j is the initial position of the i-th wolf in the j-th decision space; rand is a uniformly distributed random number in the interval [0,1].
Step 2: (Scouting) Choose the objective function f as the prey odor concentration Y t (i) and calculate the prey odor concentration at the location of the artificial wolf.According to the Y t (i), sort all artificial wolf positions in descending order, then select the first artificial wolf of the sorted population as a head wolf, whose location is regarded as x lead and the prey odor concentration is Y lead .Choose the second to S + 1 artificial wolves as the explore wolves that total S = Round(a * n), and the remaining artificial wolves as the fierce wolves.The explore wolves scout to conduct a fine search, until the maximum number of walks is reached or the maximum global optimal solution is found.
Step 3: (Summoning) The fierce wolfs quickly approach the head wolf to achieve the global convergence of solution, until the Euclidean distance d is (i) between all fierce wolves and the head wolf is less than the judgment value of distance d near or a contemporary global optimal solution is found.The calculation formula of the distance judgment value is: where d near is the judgment value of distance; D is the spatial dimension; W is the distance determining factor; max d , min d are the upper and lower limits of the decision variable x i,j , respectively.
Step 4: (Sieging) Under the command of the head wolf, other artificial wolves further update the location to achieve local fine optimization.
Step 5: (Competitive update) According to the population, update factor b, calculate the amount of the population, update R and regenerate R artificial wolves, which can be alternated with the previous generation of artificial wolves with a low fitness value to ensure the diversity of solutions and to avoid falling into the local optimal.
Step 6: Determine whether the maximum number of iterations has been reached; if it has, output the optimal global solution.If not, proceed to the next generation calculation and return to Step 2 until the maximum number of iterations is reached.
The WPA flow chart for solving the Muskingum model parameters is shown in Figure 2.
previous generation of artificial wolves with a low fitness value to ensure the diversity of solutions and to avoid falling into the local optimal.
Step 6: Determine whether the maximum number of iterations has been reached; if it has, output the optimal global solution.If not, proceed to the next generation calculation and return to Step 2 until the maximum number of iterations is reached.
The WPA flow chart for solving the Muskingum model parameters is shown in Figure 2.

Parameter Sensitivity Analysis
The WPA involves a relatively large number of parameters; the main sensitivity parameters are the distance determination factor W and the step factor S .Therefore, this paper focuses on the discussion of parameters W and S on the algorithm.In order to determine the distance determination factor W and the step factor S , the parameters were set as follow: Artificial wolf population size

Parameter Sensitivity Analysis
The WPA involves a relatively large number of parameters; the main sensitivity parameters are the distance determination factor W and the step factor S. Therefore, this paper focuses on the discussion of parameters W and S on the algorithm.In order to determine the distance determination factor W and the step factor S, the parameters were set as follow: Artificial wolf population size n = 20, explore wolf proportion factor a = 4, population update factor b = 2, the maximum number of scouting T max = 2, the number of iterations gen = 30.With the Muskingum principle, the algorithm is used to solve D-dimensional space problem (D = 2), and the upper and lower limits of the variables are max d = 1 and min d = −1, respectively.The WPA's raid step: where step d b is the step size of the d-dimensional space; S is the step factor; max d , min d are the upper and lower limits of the decision variable x i,j , respectively.
In the summoning act, the condition of the artificial wolf raid termination is d is < d near , so the rapid step should meet the following formula: where the parameters are the same as above.
From Formula (10): S > 4W (11) where the parameters are the same as above.
Water 2018, 10, 1415 7 of 14 From the above derivation, it can be seen that the normal operation of the algorithm can be satisfied when the step factor S and the distance determination factor W satisfy the formula (11).
Fixed distance determination factor W = 20, set S = 100, 105, 110, 115, 120, and 130, respectively.For each S, run the program code 20 times independently, select the evaluation index of the algorithm, the absolute deviation of the flood process, and the observed and simulated flood peak deviation.The absolute deviation of the flood process and the formula are as follow: where δ is the absolute deviation of the flood process, m 3 /s.The flood peak deviation is: where f is the flood peak deviation; Q peak obs,i and Q peak sim,i are the observed and simulated maximum outflow at peak flow event number i, respectively; n is the number of simulations times.
The evaluation index values are taken from the average of 20 independent runs.Taking the flood event of 20100821 as an example, the results are shown in Table 3.As shown in Table 3: (1) With the increase of the S value, δ decreases slightly from 7892.7 m 3 /s to 7891.6 m 3 /s, which indicates that the larger the step size, the finer the search, and the closer the flood forecast results are to the actual flood process.
(2) When W = 20, the variation range of δ and flood peak are not large, which indicates that the step size of the WPA has little influence on the forecast results in a certain range.
Keeping in mind that the smaller the searching step, the more time consuming the process is, the initial selection of S = 120, that is S = 6W.W = 20, 30, 40, 50, 60, and 70 are selected to simulate the same flood, and the results are shown in Table 4.
As shown in Table 4: (1) With the increase of W, δ decreases slightly, and the minimum value appears at W = 60.The peak deviation is also minimized, and then the W value increases; the optimization effect has no obvious change.
(2) The range of δ is not significant, and indicates that the parameter W of the wolves algorithm on the two-dimensional nonlinear optimization problem has little influence on the algorithm, and the algorithm is more applicable.With the increase of W, the raid step will be smaller, and the optimization result is too fine, which will cause the artificial wolf to turn from difficult to besieging.Thus, the algorithm has the possibility of entering an infinite loop.
In order to further analyze the influence of different weight coefficients on the forecast results, we set different weight coefficients using the same method.The results are shown in Table 5.As shown in Table 5: (1) With the decrease of φ 1 , peak deviation increased slightly, with an increase 0.1 m 3 /s only.When φ 1 is 0.8 or 0.9, the flood peak deviation takes the minimum, which indicated that the weight factor of Objective 1 focuses on the minimum deviation of flood peak.
(2) With the increase of φ 2 , the absolute deviation of the flood process is gradually reduced from 7916.9 m 3 /s to 7891.6 m 3 /s, and reaches the minimum value when φ 2 is 0.7, which indicated that the weight factor of Objective 2 focuses on the total error of the entire flood process.
According to comprehensive analysis, the weight coefficient φ 2 had a great influence on the fitting effect of the whole flood process and φ 1 had little effect on the flood peak simulation.Therefore, φ 1 = 0.3 and φ 2 = 0.7 were selected for the multi-objective weight value.

Results and Discussion
In this paper, we took flood forecasting of Ankang hydropower station to Ankang city and Ankang hydropower station to Shuhe hydropower Station as examples.Considering the length of the two Reaches, the actual demand, and the characteristics of the inflow, the calculation time steps of Reach 1 and Reach 2 were 1 h and 6 h, respectively.The improved Muskingum method (WPA), Particle Swarm Optimization (PSO), and Trial Algorithm (TA) were selected to compare to the observed outflow.In the WPA, the values of distance determination (W) and the step (S) were both 100.The maximum number of iterations was 50.The value of the artificial wolf population size (n) was 20.The values of the explore wolf proportion (a) and population update (b) were 4 and 2, respectively.In the TA [22], the parameters (k, x) of this algorithm are shown in Table 6.In the PSO algorithm, a population of 50 individuals was used.The maximum number of iterations in the program was 100.The values of the acceleration constants c 1 and c 2 were both set to 0.2.
The corresponding parameters estimated by the WPA, TA, and PSO are listed in Table 6.It can be seen from Table 6 that the parameters obtained from the three methods are all technically reasonable, but do show marked differences.The estimated value of the parameters differ very little between the WPA and PSO, but the parameters of TA vary quite significantly with WPS and PSO.The comparison of the observed and simulation results of these three methods is presented in Figures 3 and 4.
As shown in Figures 3 and 4: (1) No matter the length of the river, the simulation results of the WPA and PSO compared to those of the TA have a high fitting effect.It is embodied in the following aspects: The simulated peak values are the same as the observed ones, and the simulation process of the backwater segment is almost coincident with the observed process.However, the simulated values by the WPA and PSO in the rising water segment are different from the observed ones.Compared with the observed values, the deviation errors of the simulation values of the WPA are smaller than those of the PSO, and the deviation errors of the simulation values of the TA are the biggest.
(2) The fitting degree of the flood process in short reach is higher than that of the long river section by the WPA, PSO, and TA.Regardless of the length of the river, the error of the TA is the largest, and the WPA simulation deviation is slightly lower than that of the PSO.
In order to visualize the fitting effect of the three methods, the 20140909 flood of Ankang Power Station to Shuhe Power Station was selected as an example (Table 7), and the correlation diagram of observed and simulated flow are drawn in Figure 5.
As shown in Figure 5: (1) The flooding process of the WPA simulation has 4 points falling on the line, and the trend line coefficient of WPA simulated flood process is 0.99, which is close to 1 compared to that of the PSO and TA, which indicates that the WPA has a high precision in the simulation of the whole flood process.
(2) The simulated flood peeks of the WPA fall in the straight line y = x, which indicates that the simulation value is equal to the observed.The peak values of the PSO and TA fall on the upper side of the line y = x, which indicates that the WPA has better simulation effect on peak value.
In order to further assess the precision of the WPA on the flood process and flood peak transmission time, the absolute deviation (δ) of flood process, flood peak deviation, and flood peak transmission time error are selected as the evaluation indices.Tables 8 and 9 show the indices values by the WPA, PAO, and TA of different floods.As shown in Tables 8 and 9: (1) The WPA and PSO simulation results have a significant advantage over the results of the TA.Compared to that of the TA, the depreciation of absolute deviation cumulative value of flood process calculated by the WPA and PSO are 53% and 52%, respectively.The depreciation of absolute deviation of flood peak calculated by the WPA and PSO are 99% and 77% based on the 20140909 flood of Ankang power station-Shuhe hydropower station, respectively.As shown in Tables 8 and 9: (1) The WPA and PSO simulation results have a significant advantage over the results of the TA.Compared to that of the TA, the depreciation of absolute deviation cumulative value of flood process calculated by the WPA and PSO are 53% and 52%, respectively.The depreciation of absolute deviation of flood peak calculated by the WPA and PSO are 99% and 77% based on the 20140909 flood of Ankang power station-Shuhe hydropower station, respectively.
(2) The maximum deviations of the WPA, PSO, and TA are 9.5 m 3 /s, 197.1 m 3 /s and 865.8 m 3 /s, respectively.Moreover, the minimum deviations are 0.3 m 3 /s, 15 m 3 /s, and 402 m 3 /s, respectively, which indicate that the flood peak simulated by the WPA and PSO has obvious advantages over that of the TA.In the 20140909 flood simulation, flood peak deviations are 9.5 m 3 /s by the WPA and 197.1 m 3 /s by the PSO, respectively.Flood peak deviation simulated by the WPA is much smaller than that of the PSO, which indicates that the WPA can significantly improve peak forecast accuracy.(2) The maximum deviations of the WPA, PSO, and TA are 9.5 m 3 /s, 197.1 m 3 /s and 865.8 m 3 /s, respectively.Moreover, the minimum deviations are 0.3 m 3 /s, 15 m 3 /s, and 402 m 3 /s, respectively, which indicate that the flood peak simulated by the WPA and PSO has obvious advantages over that of the TA.In the 20140909 flood simulation, flood peak deviations are 9.5 m 3 /s by the WPA and 197.1 m 3 /s by the PSO, respectively.Flood peak deviation simulated by the WPA is much smaller than that of the PSO, which indicates that the WPA can significantly improve peak forecast accuracy.

Conclusions
In this paper, a new multi-objective Muskingum model was established and solved by a new algorithm called the WPA, and its performance was verified by using the Ankang to Shuhe section of Hanjing River datasets.The following work has been done: (1) After a brief literature review, a novel parameter estimation approach based on the WPA was proposed for parameter estimation of the nonlinear Muskingum flood, which considered two calibration objectives in the calibration procedure: (1) The weighted sum of absolute difference between the routed and observed outflows, and (2) the coefficient of determination between the routed and observed outflows.
(2) The proposed approach is compared to the other methods (TA, PSO) for an example case from the Hanjiang River, and the results demonstrate that the WPA can achieve a high degree of accuracy to estimate the parameters and results in accurate predictions of outflow.
(3) In this study, however, the parameters of WPA, such as the distance judgment factor and the step factor, have limited its application range.A possible reason is that different rivers have its applicable parameters.In future research, more floods in different rivers basins will be used to expand its range of applications.

Figure 1 .
Figure 1.Map of the Hanjiang River and the location of major reservoirs.

Figure 1 .
Figure 1.Map of the Hanjiang River and the location of major reservoirs.

Figure 2 .
Figure 2. Flow chart of parameter estimation.
principle, the algorithm is used to solve D -dimensional space problem ( 2) = D, and the upper and lower limits of the variables are max 1 = is the step size of the d-dimensional space; S is the step factor; max d , min d are the upper and lower limits of the decision variable , i j x , respectively.In the summoning act, the condition of the artificial wolf raid termination is < is near d d , so the rapid step should meet the following formula:

Figure 2 .
Figure 2. Flow chart of parameter estimation.

Figure 4 .
Figure 4. Simulation results of different floods in Reach 2.

Figure 4 .
Figure 4. Simulation results of different floods in Reach 2.

Figure 5 .
Figure 5. Correlation of observed and simulation data.

Figure 5 .
Figure 5. Correlation of observed and simulation data.

Table 1 .
Data of river section.

Table 2 .
Data of flood.

Table 1 .
Data of river section.

Table 2 .
Data of flood.

Table 3 .
Sensitivity analysis of step factor.

Table 4 .
Sensitivity analysis of distance judgment factor.

Table 6 .
Parameter estimation results obtained from different methods.

Table 7 .
20140909 flood simulation of Ankang Power Station to Shuhe Power Station.

Table 7 .
20140909 flood simulation of Ankang Power Station to Shuhe Power Station.

Table 8 .
Analysis of Ankang hydropower station-Ankang city simulation results.

Table 8 .
Analysis of Ankang hydropower station-Ankang city simulation results.