A Class of Parameter Estimation Methods for Nonlinear Muskingum Model Using Hybrid Invasive Weed Optimization Algorithm

Nonlinear Muskingum models are important tools in hydrological forecasting. In this paper, we have come up with a class of new discretization schemes including a parameter 𝜃 to approximate the nonlinear Muskingum model based on general trapezoid formulas. The accuracy of these schemes is second order, if 𝜃 ̸= 1/3 , but interestingly when 𝜃 = 1/3 , the accuracy of the presented schemegetsimprovedtothirdorder.Then,thepresentschemesaretransformedintoanunconstrainedoptimizationproblemwhichcanbesolvedbyahybridinvasiveweedoptimization(HIWO)algorithm.Finally,anumericalexampleisprovidedtoillustratetheeffectivenessofthepresentmethods.ThenumericalresultssubstantiatethefactthatthepresentedmethodshavebetterprecisioninestimatingtheparametersofnonlinearMuskingummodels.


Introduction
Flood routing plays an essential role in calculating the shape of flood waves along open river channels.In general, there are two basic methods to route floods: hydrologic routing and hydraulic routing.Hydrologic routing method is based on the storage-continuity equation whereas the hydraulic routing method is based on continuity and momentum equations.The Muskingum model, as one hydrologic routing approach, is a popular model for flood routing, whose storage depends upon the water inflow and outflow.Analysis of model features shows that classical Muskingum models can be classified into two types in form: one is linear model and the other is nonlinear model.
In the hydrological routing approach, it is very important to estimate the parameters  and  by using historical inflow and outflow records for a specific length of a reach.Generally speaking, the parameters  and  in (2) are graphically estimated by a trial-and-error procedure.If  is obtained, the values of [()+(1−)()] are calculated by using observed data in both upstream and downstream of a river and further plotted against .The particular value which generates the loop is accepted as the best estimate of .The slope of the straight line fitted through the loop derives .
In order to estimate the parameters  and , Yoon and Padmanabhan [1] suggested and discussed three methods.Chen and Yang [2] brought forward a gray-encoded accelerating genetic algorithm.Ouyang et al. [3] developed an adaptive PSO algorithm.
IWO algorithm [13] is a new bionic intelligent algorithm which simulates the spatial weed diffusion, growth, reproduction, and competitive survival of the invasive weeds [14].IWO algorithm has been widely used in a variety of optimization problems and practical engineering problems: such as multiobjective optimization problem [15], parameter estimation of chaotic systems [16], model order reduction problem [17], global numerical optimization [18], antenna arrays problem [19,20], unit commitment problem [21], optimal power flow problem [22], flow shop scheduling problem [23], traveling salesman problem [24], and economic dispatch [25].
Furthermore, we have applied the hybrid invasive weed optimization (HIWO) algorithm to the parameter estimation problems of nonlinear Muskingum model for the first time.
The remainder of the paper is organized as follows.In Section 2, we give the preliminary knowledge and the procedure of parameter derived on Muskingum models.Section 3 presents a class of new routing procedure of nonlinear Muskingum model.The basic principle of the IWO, Nelder-Mead simplex method (NMSM), and HIWO algorithms is detailed in Section 4. Section 5 shows numerical examples and analysis of results before Section 6 concludes the paper.

Muskingum Model
McCarthy in [26] analyzed the data from the Muskingum River in Ohio and identified the following relationship: where (), (), and () denote channel storage, inflow rate, and outflow rate at time , respectively;  and  denote storage time constant and weighting factor for the river.
Although the trial-and-error procedure has been used for many years, it is time consuming and entails much subjective interpretation.Therefore, in order to avoid subjective interpretations of observed data in estimating  and , a finite difference scheme is used to resulting ordinary differential equation (1), yielding where and () and () represent observed outflow discharges and inflow discharge, respectively, at time   , Δ is the time step, and  is the total time number.
As far as we know, linear Muskingum model ( 2) is not a well-fitting technique to observe data.Thus, some studies [4,27,28] proposed the following two nonlinear Muskingum models: where  is an additional parameter.Obviously, the parameter  in these nonlinear Muskingum models is no longer an approximation of the traveling time of the flood wave.
Because the nonlinear Muskingum model in ( 6) is not as popular as the one in ( 5), many authors focus on the model in ( 5) for comparison purpose.
A review of the mentioned literature in Section 1 shows that the following routing procedure has been frequently used to obtain the calculated outflow discharges Q() at time  =   (see [5,9]).
Step 1. Assume values for the three parameters, , , and .
Step 2. Calculate () using (5), where the initial outflow is the same as the initial inflow.
Step 3. Calculate the time rate of changing storage volume Δ() using the following equation: Step 4. Use the following equation to estimate the next accumulated storage: Step 5. Calculate the next outflow as follows: where Q() represents the calculated outflow discharges at time   and Ĩ( + 1) = (() + ( + 1))/2.
In a word, almost all of the above-mentioned optimization approaches used to calculate the parameters , , and  of the nonlinear Muskingum models are based on the following approximate formula: Obviously, the truncation error of this approximate formula (10) is only (Δ).Thus, it is necessary to construct a scheme of higher accuracy to approximate differential equation (1).
Based on the generalized trapezoid formula [29], this paper first develops a class of new difference schemes which contain a parameter  to approximate differential equation (1).In doing so a class of new unconstrained nonlinear optimization problems which contain a parameter  are obtained.It should be noted that the accuracy of the presented difference schemes to approximate differential equation ( 1) can be improved to third order, when  = 1/3.In other words, we can get a parameter estimation model with better accuracy, if  = 1/3.

A Class of New Routing Procedure of Nonlinear Muskingum Model
Considering that ( 5) is the most popular nonlinear Muskingum model, this paper follows the trend and focuses on the same model although other equations can be used.Firstly, we rewrite (1) as where () = () − () − ().

Hybrid Invasive Weed Optimization (HIWO) Algorithm
4.1.Invasive Weed Optimization (IWO) Algorithm.IWO algorithm is an optimization one based on swarm intelligence [13], which can steer the search process through cooperation and competition between individuals within groups.Compared with other evolutionary algorithms, IWO algorithm maintains a population-based global search strategy.Weeds refer to those strong and aggressively growing plants, which pose a serious threat to cultivated crops [31].IWO algorithm contains four core operation steps: initialization of the population, growth and reproduction, diffusion, and competitive exclusion, which will be briefly introduced in the following sections.

Population Initialization.
IWO regards for elements of algorithm design inspired by nature, which is a numerical optimization calculation method based on population.Prior to the optimization iteration, we need to set some parameters in the IWO algorithm: the initial weed population number  min , the maximum weed population number  max , the maximum seed number  max , the minimum seed number  min , the maximum generation times , the problem dimension , the nonlinear exponent , the interval step initial value  init , and the end value  final , then randomly generate  min numbers of , and ultimately calculate the adapted degree  of .
Consider the following minimized (or maximized) function optimization problems: where  min, ≤   ≤  max, ,  = 1, 2, . . ., ,  min, and  max, represent the upper bound and lower bound of the  th dimensional variable, respectively, and  represents the dimensions of the problem.Each weed in the IWO algorithm represents a solution; that is,   = [ ,1 ,  ,2 , . . .,  , ], where   denotes the  th individual in the population.Once the IWO is initialized, the algorithm will generate randomly  min numbers of individuals to form the initial population, as shown in the following formula: where,  = 1, 2, . . ., ,  = 1, 2, . . .,  min ,  represents any uniformly distributed random number between 0 and 1.

Growth and Reproduction.
In the standard IWO algorithm, each weed breeds seeds according to its fitness (survival ability).The better fitness the individual has, the more seeds it can produce.This is in line with the natural growth law; that is, there is a linear relationship between the number of seeds produced by parent generation and the adaptable degree of matrix.The worse adaptable to the environment the weeds are, the more seeds they will produce as the current research focuses on the minimum value optimization problem.The number of seeds produced by each strain of weed is decided by the following equation [14]: where   denotes the number of seeds produced per strain of weed;  denotes the fitness value of the weed;  min denotes the minimum fitness value of the weed;  max denotes the maximum fitness value of the weed;  max denotes the maximum number of seeds produced by each strain of weed;  min denotes the minimum number of seeds produced per strain of weed; floor() means the number of taking downward to obtain an integer to the nearest one.

The Spatial Diffusion.
During the optimization iteration of the IWO algorithm, the way of producing seeds is as follows: taking the parent generation as the benchmark, copy as many as   for each solution   in the parent generation, which will be added with normal random parameters,  ( ∈ [−, ]), whose mean value  = 0 and the standard deviation is , respectively, so that each strain of weed will generate   seeds.
Normal distribution probability is generally known as the Gaussian distribution one.As one of the best important probability distribution ways in mathematics, physics, engineering, and other fields, it has a significant effect on many aspects of the statistics.If the random variable  is amenable to such a probability distribution whose mathematical expectation is  and variance  2 , denoting as (,  2 ), its probability density function can be described as follows: The random variable is also named a normal random variable, with expectation  determining its position and the standard deviation  determining the amplitude of distribution.This normal distribution random parameter  is the standard normal distribution whose expectation is 0 and the variance is 1.If the specified expectation is  and the variance is , then we only need to add up  =  *  + .
In the process of iterative computing, the standard deviation of each generation changes in accordance with the following equation: where  init denotes the initial value of the interval step,  final denotes the final value of the interval step,  denotes the current value of stepsize,  denotes the maximum generation time,  denotes the current generation, and  denotes the nonlinear adjustment factor.It can be seen from formula (23) that seeds distribute around the weeds in larger step at the preliminary stage of iteration.

Competitive Exclusion.
After some generations of the algorithm evolution, the plant population will reach the maximum population  max because of rapid propagation.However, it is hoped that the number of the plants with good survival ability (fitness) is larger than that of the plants with poor survival ability (fitness).The executive way of competitive survival rule is as follows: these newborn seeds and parent weeds will be sorted in order of the fitness value after the growth and spatial distribution of each plant.After that, we have to select  max strains of plants in accordance with the fitness values from small to large (we assume that the fitness function is to compute the minimum value) and then get rid of the remaining plants.
It is quite clear that the function of the competitive exclusion operator is to ensure that the next generation of individuals is superior to the current one, as a consequence of which, the average performance of the population will improve and gradually obtain the optimal or near-optimal solutions.4.2.Nelder-Mead Simplex Method.Nelder-Mead simplex method (NMSM) [32] is also known as a derivative-free scheme method for unconstrained optimization problems.The advantages of NMSM mainly lie in its low computation complexity, fast computation speed, and strong local search ability [33].
The NMSM [32] was proposed by scholars Nelder and Mead, who tried to obtain a minimum value for solving an objective function in a multidimensional space in 1965.The NMSM has been applied to solving all kinds of optimization problems as it can obtain an accurate value for objective functions [34].Let () denote an objective function, and   denotes the group of points in the current simplex,  ∈ {1, . . ., +1}.The NMSM method is inclusive of the following steps [34].
Compute   = (  ).If   <   , accept   and terminate the iteration; otherwise, accept   and terminate the iteration.
Step 4 (contraction (1) Extract the set of reliable negative and/or positive samples   from   with the help of   .
(3) Calculate the adaptability of weeds of  0 .
(4) If the termination criterion is met, stop the iteration of the algorithm; otherwise execute the next step.
(5) Calculate the standard deviation for each generation in accordance with formula (23). (

Experimental Results and Discussion
Here, all the algorithms (BFGS, NMSM, IWO, PSO, and HIWO) are implemented by MATLAB 7.0.4 and all numerical experiments will be run on a PC with CPU Intel CORE (TM) 2 Duo T6600 2.20 GHz, with 2.00 GB of RAM and with Windows 7 as the operating system.

Experimental Results of Benchmark Functions.
In order to verify the performance of the HIWO algorithm, we have selected 20 benchmark functions as the test set.For the convenience of comparison, the functions in this paper are exactly the same as the test functions in literature [14,35].1-9 belong to multipeak functions, and each of them owns multiple local extremum, in which 7 has 100! = 9.33 + 157 local extreme points, and 10 is the simple unimodal function in 2-dimension and 3-dimension, but it is seen as the multimodal function in more dimensions.Its search process is easy to get into local optimum.These functions can test the multimodal search capability of the algorithm to test the performance of the algorithm well.Table 1 shows the characteristics of the test function, including the search area, the theoretical optimal value, and dimension.The 20 optimization functions used in the experiment are from 1 to 20 in literature [35].The five algorithms BFGS, NMSM, IWO, PSO, and LEA [35] are selected to compare performance with HIWO algorithm.The experimental data is cited from literature [35].
The number of independent running on PC for benchmark functions is 30.The function evaluations (FES) of IWO, PSO, and HIWO algorithms for functions 1-15 are 15; the FES of IWO, PSO, and HIWO algorithms for functions 16-20 are 14.The FES of the LEA algorithm for functions 1-15 and 16-20 are much more than 15 and 14, respectively.The parameter settings of the three intelligence algorithms (PSO, IWO, and HIWO) are as follows, where  is the size of population,  is the maximum number of generations, and Max is the maximum value of objective space.For the three algorithms, the  is the same,  = 5000 (from 1 to 15) or  = 500 (from 16 to 20); the setting of FES is the same.Table 2 shows the performance comparison of different algorithms for 20 functions.In Table 2, "mean" denotes the mean values of 30 independent running processes; "std" denotes the standard deviation of 30 independent running processes.It can be seen from Table 2 that BFGS algorithm is the worst in terms of calculation precision; the BFGS algorithm cannot find the optimum solution or approximate optimum solution in solving 20 functions.The calculation precision of NMSM is better than the one of BFGS, but the calculation precision of NMSM for most of the functions is not ideal, except for functions 16-19.The calculation precision of IWO is better than the one of PSO; IWO has obvious advantages in terms of calculation precision for solving 11 functions.However, PSO is better than IWO in solving 7 functions.IWO and PSO are just the same in terms of calculation precision when they solve functions 16 and 18.The LEA algorithm outperforms the aforementioned four algorithms.However, the performance of LEA is worse the one of HIWO.LEA is much greater than HIWO algorithm in the case of function evaluations; the average optimal values of the LEA algorithm are better than those of HIWO for solving 3 and 7; mean values and standard deviations of the remaining 18 functions by LEA are lagging behind those of HIWO.HIWO can find the optimum solution or approximate optimum solution in 13 functions except functions 3, 5, 6, 8, 9, 12, and 19.
Table 3 shows the significant level test ( = 0.05) paired sample Wilcoxon signed rank test result; HIWO has obvious advantages when compared the other algorithms.The statistical results of IWO, PSO, and HIWO are presented by three metrics: obvious advantages (represented by B), roughly equal (represented by E), and a distinct disadvantage (indicated by W).It is shown that the HIWO algorithm has distinct advantages in terms of calculation precision.
It can be seen from Figure 1 to Figure 20 that, compared with IWO and PSO, HIWO has faster convergence speed in solving benchmark functions.
The convergence speed of HIWO is faster than those of IWO and PSO except for 16 and 18 from Figure 1 to Figure 20.It is shown that the HIWO algorithm can improve the probability of global search and keep the good balance between global search and local search.We can conclude that HIWO is characterized by fast convergence speed, high calculation precision, and strong robustness.

Experimental Results of Muskingum Model.
In this section, we use actual observed data of flood runoff process between Chenggouwan reach and Linqing reach of the Nanyunhe river in the Haihe Basin, Tianjin, China (the length of the river segment is 83.8 km, where there is no tributary, but a levee control exists on both sides).Lifting irrigation can occur during the water delivery and flood water may discharge into the reach when rainfall is excessively high.But these situations have little effect on the flood, the routing time   duration Δ = 12 h in 1960.More detailed data can be seen in [3].Next, we will carry out numerical experiments from the following two aspects.On the one hand, for the purpose of illustrating the advantages of using the general trapezoid formula to approximate nonlinear Muskingum model (5), we use the HIWO to solve optimization problem (17) by choosing different .
It can be seen from Table 4 that the experiment results are much better, when  = 1/3, confirming the error estimation in this paper.
Obviously, given that ( 16) is a nonlinear equation about (  ), the classical Newton iteration can be applied for solving it.The numerical results are displayed in Table 5.Furthermore,        As the numerical results in Tables 5 and 6 show, HIWO performs overall the best among the five algorithms.Meanwhile, the parameters obtained by the HIWO algorithm in    1960 can be effectively applied to the changes in river outflows of flood forecast in the future.Finally, the observed and calculated flows in 1960, 1961, and 1964 are given in Figures 21, 22, and 23, respectively.Figures 21-23 show that we can calculate approximations of flood outflows, by using the HIWO algorithm to estimate the parameters of unconstrained optimization problem (17), which could be highly close to observed outflows.To sum up, the numerical results indicate that the proposed method is efficient for estimating the parameters of nonlinear Muskingum model.

Conclusions
In this paper, a class of approaches is developed for the parameter optimization of nonlinear Muskingum models.where
(i) The PSO algorithm: population size  = 20; learning factor 1 = 2 = 1.49445; the initial weight value  start = 0.9; the final weight value  end = 0.4; the maximum velocity  max = Max/5.0;the minimum velocity  min = −Max/5.0;the calculation formulation of weight factor for the current generation is  =  start − ( start −  end ) * pow(double()/, 2).(ii) The IWO algorithm: the minimum number of weed populations  min = 8; the maximum number of weed populations  max = 10; the initial step  init = 10; the final step  fin = 1−3; the minimum number of seeds  min = 2; the maximum number of seeds  max = 3; nonlinear index  = 3. (iii) The HIWO algorithm: the minimum number of weed populations  min = 4; the maximum number of weed populations  max = 5; the settings of other parameters of HIWO are the same as IWO.

Figure 21 :
Figure 21: Comparison between calculated values and observed values in 1960.

Figure 22 :Figure 23 :
Figure 22: Comparison between calculated values and observed values in 1961.
1 , V 2 , ..., V +1 .Then, go back to Step 1.4.3.The BasicSteps of HIWO Algorithm.The main steps of executing the HIWO algorithm are shown in Algorithm 2, where  is the weed populations and  is the best individual.
(22)ulate the seed number produced by each individual of parent weed   using formula(20); compute the fitness of seeds   according to formula(22)until the numbers of   and   are greater than or equal to the maximum population  max .
(7)Compute the adaptability of seeds   .(8) After sorting   and   in accordance with the fitness value, choose the number  max of the best results as the next generation weed  +1 .(9) Update the best individual   of current generation.(10) Execute NMSM algorithm and implement a precise search.(11) Let  =  + 1; turn to Step 3.

Table 1 :
The characteristics of benchmark functions.

Table 2 :
Mean value and standard deviation comparison of different algorithms for benchmark functions.
Table 4 lists the average absolute errors (AAE) and average relative errors (ARE) of calculated outflows with those of observed outflows for different  in 1960, where the AAE and the ARE are given as follows:

Table 3 :
Wilcoxon statistic result comparison by HIWO and other algorithms.

Table 4 :
Numerical results with different  for flood routing in 1960.
we have used the same parameters presented in Table5to validate the effectiveness of flood routing in 1961 and 1964, with the AAE and the ARE, respectively, listed in Table6.

Table 5 :
Result comparison of different methods in 1960.

Table 6 :
Result comparison of different methods in 1961 and 1964.