Abstract

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 , but interestingly when , the accuracy of the presented scheme gets improved to third order. Then, the present schemes are transformed into an unconstrained optimization problem which can be solved by a hybrid invasive weed optimization (HIWO) algorithm. Finally, a numerical example is provided to illustrate the effectiveness of the present methods. The numerical results substantiate the fact that the presented methods have better precision in estimating the parameters of nonlinear Muskingum models.

1. 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 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.

In the last two decades, some nonlinear optimization algorithms have been adopted to estimate the parameters , , and of the aforementioned nonlinear Muskingum model. Commonly used algorithms include genetic algorithm [4], BFGS technique [5], particle swarm optimization algorithm [6], immune clonal selection algorithm [7], parameter-setting-free harmony search algorithm [8], Nelder-Mead simplex algorithm [9], differential evolution algorithm [10], hybrid harmony search algorithm [11], and modified honey bee mating optimization algorithm [12].

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.

2. 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), yieldingwhereand 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 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 represents the calculated outflow discharges at time and .

Step 6. Repeat Steps 2–5.

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 . In other words, we can get a parameter estimation model with better accuracy, if .

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) aswhere .

Secondly, applying the generalized trapezoid formula (see [29]) to (1), we can haveCombining (12) and (11), we haveSubstituting (13) into (14), we finally obtainFurthermore, by substituting (5) into (15), we have

With the aim of estimating the parameters , , and , we can construct the following unconstrained nonlinear optimization problem:

Remark 1. By using Taylor expansion, it is easy to see that the order of above difference schemes (15) is if . Specially, for , above difference scheme (15) is of third-order accuracy (see [30]).

4. 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.

4.1.1. 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 , the maximum weed population number , the maximum seed number , the minimum seed number , the maximum generation times , the problem dimension , the nonlinear exponent , the interval step initial value , and the end value , then randomly generate numbers of , and ultimately calculate the adapted degree of .

Consider the following minimized (or maximized) function optimization problems:where , , and represent the upper bound and lower bound of the dimensional variable, respectively, and represents the dimensions of the problem.

Each weed in the IWO algorithm represents a solution; that is, , where denotes the individual in the population. Once the IWO is initialized, the algorithm will generate randomly numbers of individuals to form the initial population, as shown in the following formula:where, , , represents any uniformly distributed random number between 0 and 1.

4.1.2. 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; denotes the minimum fitness value of the weed; denotes the maximum fitness value of the weed; denotes the maximum number of seeds produced by each strain of weed; denotes the minimum number of seeds produced per strain of weed; means the number of taking downward to obtain an integer to the nearest one.

4.1.3. 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 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 , denoting as , 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 .

The method of seed production is shown as follows:where, , means the current generation. If , ; otherwise , , . denotes weed space, denotes seed space, and refers to the normal distribution random sequence.

In the process of iterative computing, the standard deviation of each generation changes in accordance with the following equation:where denotes the initial value of the interval step, 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.

4.1.4. Competitive Exclusion

After some generations of the algorithm evolution, the plant population will reach the maximum population 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 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, . The NMSM method is inclusive of the following steps [34].

Step 1 (order). The system of inequalities, , must meet the demands of the vertices.

Step 2 (reflection). Firstly, compute the reflection point of the simplex by , , where is the centroid of the best points. Then calculate . If , accept the reflected point and terminate the iteration.

Step 3 (expandsion). If , calculate the expansion point of the simplex by , ( and ). Compute . If , accept and terminate the iteration; otherwise, accept and terminate the iteration.

Step 4 (contraction). If , implement a contraction operation between and the better point of and .(a)Outside contraction: if , calculate the outside contraction point by , (). Compute . If , accept and terminate the iteration; otherwise, continue with Step 5.(b)Inside contraction: if , compute the inside contraction point by . Calculate . If , accept and end the iteration; if not, continue with Step 5.

Step 5 (shrink). Calculate the points by . And calculate , . The unordered vertices of the simplex at the next iteration consist of . Then, go back to Step 1.

4.3. The Basic Steps 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.

Algorithm 2 (HIWO algorithm). Require , , , , , , , , , , , and
Ensure (1)Extract the set of reliable negative and/or positive samples from with the help of .(2)Initialize the weed population in accordance with formula (19), let .(3)Calculate the adaptability of weeds of .(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).(6)Calculate 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 .(7)Compute the adaptability of seeds .(8)After sorting and in accordance with the fitness value, choose the number of the best results as the next generation weed .(9)Update the best individual of current generation.(10)Execute NMSM algorithm and implement a precise search.(11)Let ; turn to Step 3.

5. 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.

5.1. 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]. belong to multipeak functions, and each of them owns multiple local extremum, in which has local extreme points, and 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 to 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 are ; the FES of IWO, PSO, and HIWO algorithms for functions are . The FES of the LEA algorithm for functions and are much more than and , 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 is the maximum value of objective space. For the three algorithms, the is the same, (from to ) or (from to ); the setting of FES is the same.(i)The PSO algorithm: population size ; learning factor ; the initial weight value ; the final weight value ; the maximum velocity ; the minimum velocity ; the calculation formulation of weight factor for the current generation is .(ii)The IWO algorithm: the minimum number of weed populations ; the maximum number of weed populations ; the initial step ; the final step ; the minimum number of seeds ; the maximum number of seeds ; nonlinear index .(iii)The HIWO algorithm: the minimum number of weed populations ; the maximum number of weed populations ; the settings of other parameters of HIWO are the same as IWO.

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 . 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 and . 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 and ; 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 , , , , , , and .

Table 3 shows the significant level test () 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 and 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.

5.2. 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  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 . 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:It can be seen from Table 4 that the experiment results are much better, when , confirming the error estimation in this paper.

On the other hand, in order to test the performance of HIWO, we list the numerical results of BFGS [5], PSO [6], NMSM [9], IWO, and HIWO in 1960, where the initial parameter values of the NMSM and BFGS are chosen as and . Once the parameters , , and are calculated, we can use (16) to calculate the next calculated outflow , where . 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, we have used the same parameters presented in Table 5 to validate the effectiveness of flood routing in 1961 and 1964, with the AAE and the ARE, respectively, listed in Table 6.

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 2123 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.

6. Conclusions

In this paper, a class of approaches is developed for the parameter optimization of nonlinear Muskingum models. These methods are based on the general trapezoid formulas and a hybrid invasive weed optimization algorithm. In order to improve the precision of parameter estimation of the nonlinear Muskingum model, we first use the general trapezoid formulas to approximate the nonlinear Muskingum model and then obtain a class of high accuracy difference schemes. Last but not least, a hybrid invasive weed optimization algorithm is developed to better estimate the parameters , , and of the nonlinear Muskingum model. By virtue of the adopted Nelder-Mead simplex algorithm, the efficiency and accuracy of the new algorithm are much better compared to those of existing algorithms. The proposed algorithm has been applied to the nonlinear Muskingum model and offers similarly encouraging results.

Appendix

Considerwhere and are random integers in and is a random number in :where , , , , , , , , , , and , , , , , , , , , , where

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This work is supported by the National Science Foundation of China (11301044), the key projects of excellent young talents fund in universities of Anhui province (2013SQRL095ZD), the project supported by the Scientific Research Fund of Hunan Provincial Education Department (Grant no. 13C333), and the project supported by the Science and Technology Research Foundation of Hunan Province (Grant no. 2014GK3043).