Mathematical Problems in Engineering

Volume 2015 (2015), Article ID 573894, 15 pages

http://dx.doi.org/10.1155/2015/573894

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

^{1}College of Computer Science and Electronic Engineering, Hunan University, Changsha, Hunan 410082, China^{2}School of Information Science and Engineering, Hunan City University, Yiyang, Hunan 413000, China^{3}Department of Mathematics and Computer Science, Chizhou College, Chizhou, Anhui 247000, China^{4}College of Mathematics and Information Science, Guangxi University, Nanning, Guangxi 530004, China

Received 27 January 2015; Revised 7 May 2015; Accepted 18 May 2015

Academic Editor: Gisele Mophou

Copyright © 2015 Aijia Ouyang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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].