Application of the Elitist-Mutated PSO and an Improved GSA to Estimate Parameters of Linear and Nonlinear Muskingum Flood Routing Models

Heuristic search algorithms, which are characterized by faster convergence rates and can obtain better solutions than the traditional mathematical methods, are extensively used in engineering optimizations. In this paper, a newly developed elitist-mutated particle swarm optimization (EMPSO) technique and an improved gravitational search algorithm (IGSA) are successively applied to parameter estimation problems of Muskingum flood routing models. First, the global optimization performance of the EMPSO and IGSA are validated by nine standard benchmark functions. Then, to further analyse the applicability of the EMPSO and IGSA for various forms of Muskingum models, three typical structures are considered: the basic two-parameter linear Muskingum model (LMM), a three-parameter nonlinear Muskingum model (NLMM) and a four-parameter nonlinear Muskingum model which incorporates the lateral flow (NLMM-L). The problems are formulated as optimization procedures to minimize the sum of the squared deviations (SSQ) or the sum of the absolute deviations (SAD) between the observed and the estimated outflows. Comparative results of the selected numerical cases (Case 1–3) show that the EMPSO and IGSA not only rapidly converge but also obtain the same best optimal parameter vector in every run. The EMPSO and IGSA exhibit superior robustness and provide two efficient alternative approaches that can be confidently employed to estimate the parameters of both linear and nonlinear Muskingum models in engineering applications.


Introduction
Accurate forecasting of flood wave movement in natural river channels is extremely important for the real-time monitoring, alert and control of floods, which are effective non-engineering measures for preventing tremendous loss of lives and property. Two categories of approaches for flood routing exist: hydraulic and hydrologic methods [1]. The former routes flood by numerically solving the famous Saint-Venant equations, which usually has strict requirements for the topographical data of the investigated stream channel (such as channel cross-section and roughness) and complicated computations [2]. Conversely, the latter is based on the continuity and empirical storage equations and is more widely used in engineering applications due to its simplicity. The Muskingum flood routing model, developed by McCarthy [3], is the most frequently applied hydrologic technique.
As is known to all, the precise estimation of parameters is the key point in applying the Muskingum method for real-time flood forecasting [4]. This problem is always formulated and solved by determining the values of Muskingum parameters using historical inflow-outflow hydrograph data based on a specified optimization criterion (i.e., optimization objective). During the past decades, two types of diverse techniques have been developed to deal with the problem: traditional mathematical methods and heuristic optimization algorithms. The mathematical methods include the least-squares method (LSM) [5], the Hooke-Jeeves (HJ) pattern search in conjunction with the linear regression (HJ+LR), the conjugate gradient (HJ+CG) or the Davidon-Fletcher-Powell (HJ+DFP) algorithms [6], the nonlinear least-squares regression (NONLR) [7], the Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS) [8] and the Nelder-Mead simplex (NMS) algorithm [9]. However, most mathematical methods mentioned above inevitably have some drawbacks, such as special derivation conditions, a time-consuming quality or initial parameter assumptions. Therefore, numerous researchers focus on heuristic optimization algorithms that are characterized by fast convergence and the ability to obtain better solutions in recent decades, such as the harmony search (HS) [10], the genetic algorithm (GA) [11], the standard, improved or hybrid particle swarm optimization algorithms (PSOs) [12][13][14][15], the immune clonal selection algorithm (ICSA) [16], the differential evolution (DE) [17], and the cuckoo search (CS) algorithm [18].
The purpose of this research is to apply the newly developed elitist-mutated PSO (EMPSO) algorithm [19] and an improved gravitational search algorithm (IGSA) to solve parameter estimation problems of different forms of Muskingum models (one linear structure and two nonlinear structures). The proposed IGSA is based on the gravitational search algorithm (GSA) [20]. These two improved algorithms both have no previous applications for such issues. In the IGSA, a modified velocity updating rule and the elite strategy are introduced to enhance the global search ability and accelerate the convergence speed of the basic GSA, respectively. The experimental results of 9 widely-used standard benchmark functions with diverse properties demonstrate the global optimization abilities of the EMPSO and IGSA. The application cases verify their validity and advantages in handling parameter estimation problems of both linear and nonlinear Muskingum models.
The remainder of this paper is organized as follows: In Sect. 2, we provide the structures and the flood routing procedures of three important linear and nonlinear Muskingum models, whose structure complexities increase with the number of parameters from two to four. In Sect. 3, we briefly describe the newly developed EMPSO and the IGSA presented in this study, and then they are tested on 9 minimization benchmark functions. In Sect. 4, the EMPSO and the IGSA are successfully applied in numerical cases (three typical flood events). The results and analysis are also presented in this section. We discuss some conclusions of our research work in Sect. 5.

Muskingum Models
In previous decades, various forms of Muskingum models have been investigated [5,18,21]. Three typical linear or nonlinear Muskingum models and their corresponding flood routing equations or procedures are briefly described in this section: the original two-parameter linear Muskingum model (LMM) [3], a three-parameter nonlinear Muskingum model (NLMM) [5] and a four-parameter nonlinear Muskingum model that incorporates the lateral flow (NLMM-L) [18].

LMM
The original LMM, which is based on the basic hypothesis that the storage within a river reach is a weighted function of inflow and outflow rates, employs the following continuity and storage equations.
where S t = channel storage at time t; I t and O t = observed rates of inflow and outflow at time t, respectively; K = storage-time constant, which has a value that is similar to the flow travel time through the routing river reach; x = weighting factor, x 2 (0, 0.3] for stream channels and x 2 (0, 0.5] for reservoir storage. The finite difference solution for Eqs (1) and (2) and the flood routing procedure of LMM is given by Eqs (3)- (5).
whereÔ t = estimated outflow at time t; T = total number of time intervals; C 0 , C 1 and C 2 = three coefficients of LMM. Note that the LMM is a two-parameter (C 0 , C 1 ) Muskingum model because C 2 = 1 − C 0 − C 1 .

NLMM
However, the relationship between the channel storage S t and the weighted flow is not always and essentially linear in many river reaches; thus, the use of LMM may be inappropriate. Hence, an additional exponent parameter m was introduced to consider the effect of nonlinearity. The following form of nonlinear Muskingum model has been suggested [5].
As shown in Eq (6), the NLMM is a three-parameter (K, x and m) Muskingum model and the LMM is a particular form of NLMM with m = 1. The rate of outflow O t can be calculated by rearranging Eq (6): NLMM-L The LMM and NLMM are frequently viewed and discussed in the literature. However, they all disregard the lateral flow along the investigated reach despite the fact that lateral flow exists along many river reaches in actual flood events. Assuming that the lateral flow (Q lat ) linearly varies along a river reach and can be expressed as a ratio (α) of the inflow rate (Q lat = αI), O'Donnell [21] proposed another linear Muskingum model that consider lateral flow in 1985, it is expressed as Eqs (8) and (9).
Inspired by the above assumptions by O'Donnell [21], in 2014, Karahan and Gurarslan [18] proposed a new nonlinear Muskingum model that takes the lateral flow into consideration after the integration of continuity Eq (8) and the NLMM in Eq (6): As expressed in Eq (10), the NLMM-L is a four-parameter (K, x, m and α) Muskingum model. By rearranging Eq (10), the rate of outflow O t can be calculated using Eq (11).
Routing Procedures of the NLMM and NLMM-L In contrast to the LMM, the flood routing procedures of nonlinear Muskingum models NLMM and NLMM-L are highly complex. The routing procedures for the NLMM and the NLMM-L can be standardized using the following steps [2,8,18]: Step 1: Assume values of the Muskingum parameters (K, x and m for NLMM; K, x, m and α for NLMM-L).
Step 2: Calculate the storage amount (S t ) using Eq (6) for the NLMM and Eq (10) for the NLMM-L.
Step 3: After combining the continuity Eqs (1) and (8) with the corresponding outflow calculations in Eqs (7) and (11), the time rate of the storage change can be calculated using Eq (12) for the NLMM and Eq (13) for the NLMM-L.
Step 4: Calculate the next storage using Eq (14), where Δt is assumed to represent the unit time.
Step 5: Calculate the next estimated outflow ðÔ tþ1 Þ using Eq (15) for the NLMM and Eq (16) for the NLMM-L.Ô tþ1 ¼ Note that Eqs (15) and (16) use the observed inflow at the previous time-point (I t ) instead of the observed inflow at the current time (I t+1 ) compared with Eqs (7) and (11) because Eq (15) occasionally provides better estimated outflow ðÔ tþ1 Þ as suggested and reported by [6,8].
Step 6: Repeat Steps 2-5 for all time steps.

Two Improved Heuristic Algorithms
Elitist-mutated PSO Standard PSO and Its Developments. The PSO algorithm, originally introduced by Kennedy and Eberhart [22], is a population-based stochastic search technique inspired by the social behavior of fish schooling or bird flocking. In PSO, each individual within the swarm is called as a particle and represents a candidate solution to the optimization problem. For a D-dimensional search space, assume that X i = (x i1 , x i2 , . . ., x iD ) and V i = (v i1 , v i2 , . . ., v iD ) denote the position vector and the velocity vector of the ith particle, respectively. The best previously visited position of the ith particle and the current global best position in the swarm are recorded as pbest i = (p i1 , p i2 , . . ., p iD ) and gbest = (p g1 , p g2 , . . ., p gD ), respectively, where g is the index of the best particle in the swarm. During iterations, the swarm is manipulated according to the updating rules written as Eqs (17) and (18). Note that such process involves individual intelligence, i.e., the particles learn through their own experience (local search) and the experience of their peers (global search).
where d = 1, . . ., D represents the index for the decision variables; i = 1, . . ., pop and pop = the number of particles in the swarm; n = iteration number; r 1 and r 2 = uniformly generated random numbers in [0, 1]; c 1 and c 2 = cognitive and social parameters, respectively, which are referred to as acceleration constants. In addition, the value of velocity v n id in each iteration should be limited Eqs (17) and (18) yield the standard PSO that may have shortcomings of premature convergence and poor control of its search capability. To overcome these drawbacks, extended studies and developments were reported [23][24][25][26]. Among which, two big variations focus on modifying the model coefficients are as follows.
In 1988, Shi and Eberhart [26] introduced a linearly decreasing inertia weight parameter w n into Eq (17) to balance the local search and the global search abilities, which are expressed as Eqs (19) and (20) where w max and w min = the initial and the final inertia weights, respectively; n = current iteration number; and N = maximum iterations.
In 1999, Clerc [23] introduced the constriction factor χ to control the changes in velocity and assure better convergence of the PSO, as shown in Eqs (21) and (22).
Latest Velocity Updating Rule and Elitist Mutation Operator. A newly developed PSO, namely, elitist-mutated PSO (EMPSO), was proposed by Nagesh Kumar and Janga Reddy [19] for solving water resource problems. Two main improvements in the EMPSO algorithm are: (1) each particle calculates its velocity using the latest updating rule as shown in Eq (23), where χ is the constriction coefficient and w is the inertia weight; and (2) a new strategic mechanism called elitist mutation operator is introduced to enhance the diversity of the swarm and explore new regions in the whole search space. Pseudo-code of the EMPSO algorithm is presented in As illustrated in Fig 1a, in each iteration, the elitist mutation operator is performed on a predefined number (NM) of the worst fitness particles in the swarm. This process of random perturbation is described as follows: first, all particles are sorted in ascending order based on their fitness function and the index numbers for the respective particles are obtained (i.e. index of the sorted swarm is recorded in ASF[i], i = 1, . . ., pop); second, the elitist mutation is performed on the front NM worst particles (selected number of least ranked particles to be elitistmutated) and the respective particle position vectors are replaced with the new mutated position vectors obtained after performing variable-wise mutation on the global best position vector (p em is the mutation probability), whereas the velocity vectors of these particles are unvaried.

Improved GSA
Basic GSA. Based on the law of gravity and mass interactions, in 2009, Rashedi and Nezamabadi-pour [20] proposed a novel heuristic algorithm, namely, the gravitational search algorithm (GSA). For an optimization problem, the searcher agents in GSA are a collection of masses in which the values of the masses are proportional to their fitness functions. During the iterative process, the masses interact with each other according to Newtonian gravity and the laws of motion. A heavier mass has a higher attraction, which indicates greater efficiency (similar to the global optimum) and a slower speed of movement. The basic GSA is mathematically described as follows.
Consider a system with pop agents (masses), in which the position of the ith agent (candidate solution) is denoted by where x d i represents the position of the ith agent in the dth dimension and D is the space dimension. According to the Newton law of gravity, the force acting on the ith mass from the jth mass at time t is defined as Eq (25). The total force acting on agent i in the dimension d is considered to be a randomly weighted sum of the forces exerted from other agents, as expressed in Eq (26) where M pi , M aj = the passive and active gravitational masses of agent i and agent j, respectively; M ai = M pi = M i i = 1, . . ., pop, M i is the inertial mass of the ith agent; G(t) = gravitational constant at time t; ε = a small constant; R ij (t) = ||X i (t), X j (t)|| 2 , is the Euclidian distance between agents i and j; rand j = a uniformly generated random number in [0, 1]; Kbest = the set of first K agents with the best fitness values and the largest masses, which is a function of time t, has the initial value K 0 = pop and is linearly decreases to 1 at the end of each iteration. Note that the gravitational constant G is a function of the initial value (G 0 , a problemdependent parameter) and time t: Based on the law of motion, the acceleration of the ith agent in the dth dimension at time t is given by The next velocity of each agent i is considered to be a fraction of its current velocity added to its acceleration and is expressed as follows: where fit i (t) = the fitness value of the agent i at time t; worst(t) and best(t) = the worst fitness and the best fitness, respectively, among all agents. Modified Velocity Updating Rule and Elite Strategy. The basic GSA may spend a significant amount of time converging to the global optimum due to the presence of heavier masses at the end of every run. Therefore, we propose an improved GSA (IGSA) which employs the following two strategies to overcome this drawback. The first strategy learns from the idea of memory and social information of PSO and defines a new velocity updating rule for agents, which is written as Eq (32) [31]. The second strategy adds the elite strategy to GSA to accelerate its convergence speed. The idea is to directly preserve a certain number of elite agents in the current generation and replace an equal number of worst agents of the new generated offspring generation. The top 5% of agents are preserved in each generation. Pseudo-code of the IGSA is shown in Fig 2. v where rand i , r 1 and r 2 = uniformly generated random numbers in [0, 1]; c 1 and c 2 = weighting factors; and x d g ðtÞ = the current best solution.

Performance Test Using Benchmark Functions
Benchmark functions are commonly recognized as an important tool to validate the performance of optimization algorithms [24,25,32]. There have been many kinds of benchmark functions reported in the literature [33,34]. However, only a comprehensive selection of benchmark functions with various characteristics can be truly useful to test new algorithms in an unbiased way. For this reason, a rich test suite of 9 standard minimization benchmark functions with diverse properties in terms of separability, modality were used as experiments for evaluating the EMPSO and IGSA. Table 1 gives a detailed description of these functions, where D is the dimension of the function, f min is the optimum value of the function. Functions f1-f7 are high-dimensional problems. The first four functions (f1 to f4) are unimodal, The values of a ij in f8 are given in S1 Table. The vectors a i and c i in f9 are given in S2 Table. doi:10.1371/journal.pone.0147338.t001 which are relatively easy to solve. Functions f5-f9 are multimodal so that the algorithm really suffers from being premature. Functions f5-f7, where the number of local minima increases exponentially with the problem dimension, appear to be the most difficult class of optimization problems. Functions f8 and f9 are two low-dimensional functions which have only a few local minima. We applied the EMPSO and IGSA to the above 9 benchmark functions and compared the experimental results with those obtained by the standard PSO as well as basic GSA. The results are averaged over 50 independent runs and the best-so-far solution, mean and standard deviation of the best solution in each run are listed in Table 2  Best: best-so-far solution over 50 runs.
Mean: mean of the best solutions in 50 runs. Std.: standard deviation of the best solutions in 50 runs.
As can be seen from Table 2, the best results are indicated in bold font. Generally speaking, the EMPSO provides much better results than PSO for all the 9 benchmark functions according to the three statistics (Best, Mean and Std.). The IGSA can find better solutions than GSA for functions f1-f7 and it strikingly improves the robustness of GSA (smaller values of Std.) on all functions except for function f6. If comprehensively consider the values of Best and Std., the EMPSO performs the best on 5 functions (f1, f2, f4, f6, f9) and IGSA is the best on the other 4 functions (f3, f5, f7, f8). The results in Table 2 also show that EMPSO and IGSA have better global optimization abilities than the PSO and GSA in solving most of the 9 benchmark functions and can obtain similar solutions.

Numerical Cases
In the parameter optimization problems of Muskingum models, minimization of the sum of the squared deviations (SSQ) or the sum of the absolute deviations (SAD) between the observed and the estimated outflows is always adopted as the objective function f, defined as follows: where O t = observed outflow at time t;Ô t ðPÞ = estimated outflow at time t by the Muskingum routing equation that is Eq (4) for the LMM, Eq (15) for the NLMM and Eq (16) for the NLMM-L; P = the parameter vector need to be calibrated, where P = (C 0 , C 1 ) in LMM, P = (K, x, m) in NLMM, P = (K, x, m, α) in NLMM-L.
To evaluate the practicability of the EMPSO algorithm and the IGSA in engineering applications, we applied these two improved heuristic algorithms to seek the optimal parameter vector P for the three different Muskingum models and compared the results with those obtained by RGA and standard PSO, as well as the basic GSA. The optimal parameter vectors obtained in this study are also compared with the best existing solutions reported in previous literature. For the above five algorithms, the iterations proceed until the stopping criterion is satisfied, which is expressed as where n is the iteration number and N is the maximum number of iterations (set to 5000); f best (n) is the best value of f in the nth iteration and δ is convergence accuracy.
For the five algorithms in applications, the population size pop was set to 50 and they were implemented on a PC with a 32-bit Windows 7 operating system, 4 GB RAM and 2.93 GHzcore (TM) i3-based processor. Each algorithm was performed over 50 runs on the three Muskingum models for the numerical examples. In RGA, the tournament selection, simulated binary crossover (SBX) and polynomial mutation operators [35] are used; the crossover probability p c = 0.85 and mutation probability p m = 0.05; the distribution index for SBX is 10 and the distribution index for the mutation operator is 100. Parameter settings of the PSO, EMPSO, GSA and IGSA are the same with Sect. 3.

Case 1: Application to LMM
Flood Data from the South Canal of China in August 1961. A flood occurred in the south canal of China between the Linqing River and the Chenggou Bay in August 1961, in which the inflow and outflow hydrographs exhibit obvious linear characteristics [36], is employed as the numerical case for the LMM, where Δt = 12h and T = 28. The search ranges for the two parameters in LMM are set to C 0 , C 1 2 (0.00, 0.50). One best existing solution according to the literature [11] is C 0 = 0.4736, C 1 = 0.0301 and SAD = 141.225, which is used as a reference.
Results and Analysis. For comparison, the statistical results (Best, Worst, Mean, and Std.) of the SAD, the model parameters, the iteration number and the CPU time for convergence (convergence accuracy δ) by the five algorithms (RGA, PSO, GSA, EMPSO, IGSA) are listed in Table 3: (1) with the exception of RGA, the other four algorithms find the same optimal solution (SAD = 141.194; C 0 = 0.4729 and C 0 = 0.0317) after 50 runs, which is better than the reference; however, only the GSA, EMPSO and IGSA can steadily converge to the same optimal solution for the LMM in every run (values of Std. for the SADs and the parameters are 0.0000E +00), whereas the optimal solutions obtained by the PSO are slightly fluctuate between different runs; (2) the GSA and IGSA require more time to converge than other three algorithms; (3) compared with GSA and IGSA, the EMPSO has a faster convergence speed (only requires 120 iterations and an average 0.0067s of CPU time for convergence in every run) and exhibits better stability (smallest values of Std.). The estimated outflow hydrograph by the LMM using the best parameter vector obtained in this study is shown in Fig 3. Fig 4 shows the comparison of the average convergence rate among the five algorithms on the LMM.  Wilson (1974). The data set from ref. [37], which had been demonstrated to have a nonlinear relationship between the storage and the weighted-flow [7], is taken as the numerical case for the NLMM. It is a single peak hydrograph that has been previously investigated by many researchers [2, 4, 7, 8, 10, 16-18, 21, 38], where Δt = 6h and T = 21. The NLMM has three parameters and their search ranges are set to K 2 [0.01, 1.00], x 2 [0.00, 0.30]and m 2 [1.00, 3.00]. One best existing solution for the NLMM refers to Xu and Qiu [17] using the differential evolution (DE) algorithm was K = 0.5175, x = 0.2869, m = 1.8680 and SSQ = 36.77.
Results and Analysis. The statistical results, which resemble Table 3, are also listed in Table 4: (1) only the EMPSO and IGSA can steadily converge to the same optimal solution (SSQ = 36.7679; K = 0.5175, x = 0.2869 and m = 1.8680) for the NLMM in every run, whereas optimal solutions obtained by RGA, PSO and GSA fluctuate between different runs; (2) for GSA, EMPSO and IGSA, the average number of iterations that are required for convergence in every run are approximately 2062, 191 and 780, respectively, and the corresponding CPU consuming times are 0.9992s, 0.0883s and 0.4864s; these data indicate that the EMPSO has a faster convergence speed and the IGSA obviously improves the convergence performance of GSA; (3) the EMPSO has the best performance in optimizing the NLMM (the lowest CPU time for every run and a steady fluctuation among of iterations with the lowest Std. = 51.22). The estimated outflow hydrograph by the NLMM using the best parameter vector obtained in this   Table 4. Statistics of different algorithms performed on the NLMM over 50 runs for the data set of Wilson (1974).  Results and Analysis. For this numerical case to the NLMM-L, statistical results resemble Table 3 are presented in Table 5. As shown in Table 5: (1) only the EMPSO and IGSA can steadily find the global optimal solution in every run (the values of Std. for the SSQs and the parameters are 0.0000E+00), which is same to the reference; (2) the EMPSO still has the best performance in optimizing the NLMM-L (the lowest average CPU time for every run and a fairly steady fluctuation among iterations with the lowest Std. = 39.82). The estimated outflow hydrograph by the NLMM-L using the best solution obtained in this study is shown in Fig 7. Fig 8 shows the average convergence rate between the five different algorithms on the NLMM-L.

Conclusions
In this study, the EMPSO algorithm and the IGSA were applied for solving the parameter estimation problems of three forms of linear or nonlinear Muskingum models (LMM, NLMM and NLMM-L). The LMM has two parameters and the NLMM has three, whereas the NLMM-L considers the lateral flow along the river reach, which has a more complex structure with four parameters. The EMPSO and IGSA were tested on a rich set of 9 standard minimization benchmark functions. Then three typical flood events used in previous literature were selected as numerical cases (Case 1-3) to evaluate the practicability of the EMPSO and IGSA in applications. The results by the EMPSO and IGSA were compared with those obtained by the RGA, PSO and GSA, as well as the best reported solutions in the literature. Several conclusions are summarized as follows.
1. only the EMPSO and IGSA can steadily converge to the same optimal solution for the three Muskingum models in every run compared with RGA, standard PSO and the basic GSA; 2. the GSA may require more iterations and CPU time than the IGSA to find the same optimal solution for the LMM, and the results obtained by the GSA for the NLMM and the NLMM-L are the worst among the five algorithms (the largest values of SSQ and Std. in Tables 4 and 5), which indicates that the proposed IGSA can improve the performance (including the search efficiency, convergence speed and the stability) of GSA in optimizing the NLMM and the NLMM-L; 3. the EMPSO has the fastest convergence rate and the best robustness than the other four algorithms for the three Muskingum models in term of specified convergence accuracy and the Std. values.
Supporting Information S1