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

Authors: Aijia Ouyang; Li-Bin Liu; Zhou Sheng; Fan Wu
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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 θ
≠
1
/
3, but interestingly when θ
=
1
/
3, 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.


---

## Body

## 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 parametersK and x by using historical inflow and outflow records for a specific length of a reach. Generally speaking, the parameters K and x in (2) are graphically estimated by a trial-and-error procedure. If x is obtained, the values of [
x
I
(
t
)
+
(
1
-
x
)
Q
(
t
)
] are calculated by using observed data in both upstream and downstream of a river and further plotted against S. The particular value which generates the loop is accepted as the best estimate of x. The slope of the straight line fitted through the loop derives K.In order to estimate the parametersK and x, 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 parametersK, x, and m 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 Section2, 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:(1)
d
S
t
d
t
=
I
t
-
Q
t
,
(2)
S
t
=
K
x
I
t
+
1
-
x
Q
t
,where S
(
t
), I
(
t
), and Q
(
t
) denote channel storage, inflow rate, and outflow rate at time t, respectively; K and x 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 estimatingK and x, a finite difference scheme is used to resulting ordinary differential equation (1), yielding(3)
S
i
-
S
i
-
1
Δ
t
=
I
i
+
I
i
-
1
2
-
Q
i
+
Q
i
-
1
2
,where(4)
S
i
=
K
x
I
i
+
1
-
x
Q
i
,
S
i
-
1
=
K
x
I
i
-
1
+
1
-
x
Q
i
-
1
,
i
=
2
,
…
,
n
,and Q
(
i
) and I
(
i
) represent observed outflow discharges and inflow discharge, respectively, at time t
i, Δ
t is the time step, and n 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:(5)
S
t
=
K
x
I
t
+
1
-
x
Q
t
m
,
(6)
S
t
=
K
x
I
m
t
+
1
-
x
Q
m
t
,where m is an additional parameter. Obviously, the parameter K 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 Section1 shows that the following routing procedure has been frequently used to obtain the calculated outflow discharges Q
~
(
i
) at time t
=
t
i (see [5, 9]).Step 1.
Assume values for the three parameters,K, x, and m.Step 2.
CalculateS
(
i
) using (5), where the initial outflow is the same as the initial inflow.Step 3.
Calculate the time rate of changing storage volumeΔ
S
(
i
) using the following equation:(7)
Δ
S
i
Δ
t
=
-
1
1
-
x
S
i
K
1
/
m
+
1
1
-
x
I
i
.Step 4.
Use the following equation to estimate the next accumulated storage:(8)
S
i
+
1
=
S
i
+
Δ
S
i
.Step 5.
Calculate the next outflow as follows:(9)
Q
~
i
+
1
=
1
1
-
x
S
i
+
1
K
1
/
m
-
x
1
-
x
I
~
i
+
1
,where Q
~
(
i
) represents the calculated outflow discharges at time t
i and I
~
(
i
+
1
)
=
(
I
(
i
)
+
I
(
i
+
1
)
)
/
2.Step 6.
Repeat Steps 2–5.In a word, almost all of the above-mentioned optimization approaches used to calculate the parametersK, x, and m of the nonlinear Muskingum models are based on the following approximate formula:(10)
Δ
S
t
i
Δ
t
≈
d
S
t
i
d
t
=
-
1
1
-
x
S
t
i
K
1
/
m
+
1
1
-
x
I
t
i
.Obviously, the truncation error of this approximate formula (10) is only O
Δ
t. 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.


## 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(11)
d
S
t
d
t
=
S
t
+
p
t
,where p
(
t
)
=
I
(
t
)
-
Q
(
t
)
-
S
(
t
).Secondly, applying the generalized trapezoid formula (see [29]) to (1), we can have(12)
S
~
i
=
S
i
+
1
-
Δ
t
d
S
i
+
1
d
t
,
S
i
+
1
=
S
i
+
Δ
t
2
1
-
θ
d
S
i
d
t
+
θ
d
S
~
i
d
t
+
d
S
i
+
1
d
t
.Combining (12) and (11), we have(13)
S
~
i
=
1
-
Δ
t
S
i
+
1
-
Δ
t
p
i
+
1
,
(14)
S
i
+
1
=
S
i
+
Δ
t
2
p
i
+
1
+
p
i
+
1
-
θ
S
i
+
θ
S
~
i
+
S
i
+
1
.Substituting (13) into (14), we finally obtain(15)
1
-
Δ
t
2
θ
S
i
+
1
=
1
-
Δ
t
2
θ
S
i
+
Δ
t
2
I
t
j
-
Q
i
+
Δ
t
2
1
-
Δ
t
θ
I
i
+
1
-
Q
i
+
1
.Furthermore, by substituting (5) into (15), we have(16)
1
-
Δ
t
2
θ
K
x
I
i
+
1
+
1
-
x
Q
i
+
1
m
=
1
-
Δ
t
2
θ
K
x
I
i
+
1
-
x
Q
i
m
+
Δ
t
2
I
i
-
Q
i
+
Δ
t
2
1
-
Δ
t
θ
I
i
+
1
-
Q
i
+
1
.With the aim of estimating the parametersK, x, and m, we can construct the following unconstrained nonlinear optimization problem:(17)
m
i
n
f
K
,
x
,
m
=
∑
i
=
1
n
-
1
1
-
Δ
t
2
θ
K
x
I
i
+
1
+
1
-
x
Q
i
+
1
m
-
1
-
Δ
t
2
θ
K
x
I
i
+
1
-
x
Q
i
m
-
Δ
t
2
I
i
-
Q
i
+
Δ
t
2
1
-
Δ
t
θ
I
i
+
1
-
Q
i
+
1
.Remark 1.
By using Taylor expansion, it is easy to see that the order of above difference schemes (15) is O
(
Δ
t
)
2 if θ
≠
1
/
3. Specially, for θ
=
1
/
3, 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 numberP
m
i
n, the maximum weed population number P
m
a
x, the maximum seed number s
m
a
x, the minimum seed number s
m
i
n, the maximum generation times G, the problem dimension d, the nonlinear exponent n, the interval step initial value σ
i
n
i
t, and the end value σ
f
i
n
a
l, then randomly generate P
m
i
n numbers of W, and ultimately calculate the adapted degree f of W.Consider the following minimized (or maximized) function optimization problems:(18)
m
i
n
m
a
x
f
x
,
x
=
x
1
,
x
2
,
…
,
x
d
,where x
m
i
n
,
i
≤
x
i
≤
x
m
a
x
,
i, i
=
1,2
,
…
,
d, x
m
i
n
,
i and x
m
a
x
,
i represent the upper bound and lower bound of the i
t
h dimensional variable, respectively, and d represents the dimensions of the problem.Each weed in the IWO algorithm represents a solution; that is,x
i
=
[
x
i
,
1
,
x
i
,
2
,
…
,
x
i
,
d
], where x
i denotes the i
t
h individual in the population. Once the IWO is initialized, the algorithm will generate randomly P
m
i
n numbers of individuals to form the initial population, as shown in the following formula:(19)
x
i
,
k
0
=
x
m
i
n
,
k
+
r
·
x
m
a
x
,
k
-
x
m
i
n
,
k
,where, k
=
1,2
,
…
,
d, i
=
1,2
,
…
,
P
m
i
n, r 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]:(20)
N
s
=
f
l
o
o
r
f
m
a
x
-
f
f
m
a
x
-
f
m
i
n
s
m
a
x
-
s
m
i
n
+
s
m
i
n
,where N
s denotes the number of seeds produced per strain of weed; f denotes the fitness value of the weed; f
m
i
n denotes the minimum fitness value of the weed; f
m
a
x denotes the maximum fitness value of the weed; s
m
a
x denotes the maximum number of seeds produced by each strain of weed; s
m
i
n denotes the minimum number of seeds produced per strain of weed; f
l
o
o
r
(
x
) 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 asN
s for each solution x
i in the parent generation, which will be added with normal random parameters, r
(
r
∈
[
-
σ
,
σ
]
), whose mean value μ
=
0 and the standard deviation is σ, respectively, so that each strain of weed will generate N
s 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 variabler is amenable to such a probability distribution whose mathematical expectation is μ and variance σ
2, denoting as N
(
μ
,
σ
2
), its probability density function can be described as follows:(21)
f
x
=
1
σ
2
π
e
x
p
-
x
-
μ
2
2
σ
2
.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 r 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 r
=
r
∗
σ
+
μ.The method of seed production is shown as follows:(22)
S
i
,
j
,
k
g
=
W
i
,
k
g
+
R
i
,
j
,
k
g
,where, i
=
1,2
,
…
,
P, g means the current generation. If g
=
1, P
=
P
m
i
n; otherwise P
=
P
m
a
x, j
=
1,2
,
…
,
N
s, k
=
1,2
,
…
,
d. W denotes weed space, S denotes seed space, and R 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:(23)
σ
=
G
-
g
G
n
σ
i
n
i
t
-
σ
f
i
n
a
l
+
σ
f
i
a
n
l
,where σ
i
n
i
t denotes the initial value of the interval step, σ
f
i
n
a
l denotes the final value of the interval step, σ denotes the current value of stepsize, G denotes the maximum generation time, g denotes the current generation, and n 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 populationP
m
a
x 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 P
m
a
x 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 f
(
x
) denote an objective function, and x
i denotes the group of points in the current simplex, i
∈
{
1
,
…
,
N
+
1
}. The NMSM method is inclusive of the following steps [34].Step 1 (order).
The system of inequalities,f
(
x
1
)
≤
f
(
x
2
)
≤
⋯
≤
f
(
x
i
)
≤
⋯
≤
f
(
x
N
+
1
), must meet the demands of the N
+
1 vertices.Step 2 (reflection).
Firstly, compute the reflection pointx
R of the simplex by x
R
=
x
¯
+
α
(
x
¯
-
x
N
+
1
), (
α
>
0
), where x
¯
=
(
1
/
N
)
∑
i
=
1
N is the centroid of the n best points. Then calculate f
R
=
f
(
x
R
). If f
1
≤
f
R
≤
f
N, accept the reflected point x
R and terminate the iteration.Step 3 (expandsion).
Iff
R
<
f
1, calculate the expansion point x
E of the simplex by x
E
=
x
¯
+
β
(
x
R
-
x
¯
), (β
>
1 and β
>
α). Compute f
E
=
f
(
x
E
). If f
E
<
f
R, accept x
E and terminate the iteration; otherwise, accept x
R and terminate the iteration.Step 4 (contraction).
Iff
R
≥
f
N, implement a contraction operation between x
¯ and the better point of x
R and x
n
+
1.(a)
Outside contraction: iff
N
≤
f
R
<
f
n
+
1, calculate the outside contraction point x
C by x
C
=
x
¯
+
γ
(
x
R
-
x
¯
), (0
<
γ
<
1). Compute f
C
=
f
(
x
C
). If f
C
<
f
R, accept x
C and terminate the iteration; otherwise, continue with Step 5.
(b)
Inside contraction: iff
R
≥
f
N
+
1, compute the inside contraction point x
C
C by x
C
C
=
x
¯
-
γ
(
x
¯
-
x
N
+
1
). Calculate f
C
C
=
f
(
x
C
C
). If f
C
C
<
f
R, accept x
C
C and end the iteration; if not, continue with Step 5.Step 5 (shrink).
Calculate theN points by v
i
=
x
1
+
σ
(
x
i
-
x
1
)
(
0
<
σ
<
1
). And calculate f
(
v
i
), i
=
2
,
…
,
N
+
1. The unordered vertices of the simplex at the next iteration consist of x
1
,
v
2
,
…
,
v
N
+
1. 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 Algorithm2, where W is the weed populations and b is the best individual.Algorithm 2 (HIWO algorithm).
RequireW, G, P
m
a
x, P
m
i
n, d, n, σ
i
n
i
t, σ
f
i
n
a
l, s
m
a
x, s
m
i
n, x
m
i
n, and x
m
a
x
.
Ensureb
.
(1)
Extract the set of reliable negative and/or positive samplesT
n from U
n with the help of P
n.
(2)
Initialize the weed populationW
0 in accordance with formula (19), let g
=
0.
(3)
Calculate the adaptability of weeds ofW
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).
(6)
Calculate the seed number produced by each individual of parent weedW
g using formula (20); compute the fitness of seeds S
g according to formula (22) until the numbers of W
g and S
g are greater than or equal to the maximum population P
m
a
x.
(7)
Compute the adaptability of seedsS
g.
(8)
After sortingW
g and S
g in accordance with the fitness value, choose the number P
m
a
x of the best results as the next generation weed W
g
+
1.
(9)
Update the best individualb
g of current generation.
(10)
Execute NMSM algorithm and implement a precise search.
(11)
Letg
=
g
+
1; turn to Step 3.


## 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 numberP
m
i
n, the maximum weed population number P
m
a
x, the maximum seed number s
m
a
x, the minimum seed number s
m
i
n, the maximum generation times G, the problem dimension d, the nonlinear exponent n, the interval step initial value σ
i
n
i
t, and the end value σ
f
i
n
a
l, then randomly generate P
m
i
n numbers of W, and ultimately calculate the adapted degree f of W.Consider the following minimized (or maximized) function optimization problems:(18)
m
i
n
m
a
x
f
x
,
x
=
x
1
,
x
2
,
…
,
x
d
,where x
m
i
n
,
i
≤
x
i
≤
x
m
a
x
,
i, i
=
1,2
,
…
,
d, x
m
i
n
,
i and x
m
a
x
,
i represent the upper bound and lower bound of the i
t
h dimensional variable, respectively, and d represents the dimensions of the problem.Each weed in the IWO algorithm represents a solution; that is,x
i
=
[
x
i
,
1
,
x
i
,
2
,
…
,
x
i
,
d
], where x
i denotes the i
t
h individual in the population. Once the IWO is initialized, the algorithm will generate randomly P
m
i
n numbers of individuals to form the initial population, as shown in the following formula:(19)
x
i
,
k
0
=
x
m
i
n
,
k
+
r
·
x
m
a
x
,
k
-
x
m
i
n
,
k
,where, k
=
1,2
,
…
,
d, i
=
1,2
,
…
,
P
m
i
n, r 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]:(20)
N
s
=
f
l
o
o
r
f
m
a
x
-
f
f
m
a
x
-
f
m
i
n
s
m
a
x
-
s
m
i
n
+
s
m
i
n
,where N
s denotes the number of seeds produced per strain of weed; f denotes the fitness value of the weed; f
m
i
n denotes the minimum fitness value of the weed; f
m
a
x denotes the maximum fitness value of the weed; s
m
a
x denotes the maximum number of seeds produced by each strain of weed; s
m
i
n denotes the minimum number of seeds produced per strain of weed; f
l
o
o
r
(
x
) 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 asN
s for each solution x
i in the parent generation, which will be added with normal random parameters, r
(
r
∈
[
-
σ
,
σ
]
), whose mean value μ
=
0 and the standard deviation is σ, respectively, so that each strain of weed will generate N
s 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 variabler is amenable to such a probability distribution whose mathematical expectation is μ and variance σ
2, denoting as N
(
μ
,
σ
2
), its probability density function can be described as follows:(21)
f
x
=
1
σ
2
π
e
x
p
-
x
-
μ
2
2
σ
2
.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 r 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 r
=
r
∗
σ
+
μ.The method of seed production is shown as follows:(22)
S
i
,
j
,
k
g
=
W
i
,
k
g
+
R
i
,
j
,
k
g
,where, i
=
1,2
,
…
,
P, g means the current generation. If g
=
1, P
=
P
m
i
n; otherwise P
=
P
m
a
x, j
=
1,2
,
…
,
N
s, k
=
1,2
,
…
,
d. W denotes weed space, S denotes seed space, and R 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:(23)
σ
=
G
-
g
G
n
σ
i
n
i
t
-
σ
f
i
n
a
l
+
σ
f
i
a
n
l
,where σ
i
n
i
t denotes the initial value of the interval step, σ
f
i
n
a
l denotes the final value of the interval step, σ denotes the current value of stepsize, G denotes the maximum generation time, g denotes the current generation, and n 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 populationP
m
a
x 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 P
m
a
x 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.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 numberP
m
i
n, the maximum weed population number P
m
a
x, the maximum seed number s
m
a
x, the minimum seed number s
m
i
n, the maximum generation times G, the problem dimension d, the nonlinear exponent n, the interval step initial value σ
i
n
i
t, and the end value σ
f
i
n
a
l, then randomly generate P
m
i
n numbers of W, and ultimately calculate the adapted degree f of W.Consider the following minimized (or maximized) function optimization problems:(18)
m
i
n
m
a
x
f
x
,
x
=
x
1
,
x
2
,
…
,
x
d
,where x
m
i
n
,
i
≤
x
i
≤
x
m
a
x
,
i, i
=
1,2
,
…
,
d, x
m
i
n
,
i and x
m
a
x
,
i represent the upper bound and lower bound of the i
t
h dimensional variable, respectively, and d represents the dimensions of the problem.Each weed in the IWO algorithm represents a solution; that is,x
i
=
[
x
i
,
1
,
x
i
,
2
,
…
,
x
i
,
d
], where x
i denotes the i
t
h individual in the population. Once the IWO is initialized, the algorithm will generate randomly P
m
i
n numbers of individuals to form the initial population, as shown in the following formula:(19)
x
i
,
k
0
=
x
m
i
n
,
k
+
r
·
x
m
a
x
,
k
-
x
m
i
n
,
k
,where, k
=
1,2
,
…
,
d, i
=
1,2
,
…
,
P
m
i
n, r 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]:(20)
N
s
=
f
l
o
o
r
f
m
a
x
-
f
f
m
a
x
-
f
m
i
n
s
m
a
x
-
s
m
i
n
+
s
m
i
n
,where N
s denotes the number of seeds produced per strain of weed; f denotes the fitness value of the weed; f
m
i
n denotes the minimum fitness value of the weed; f
m
a
x denotes the maximum fitness value of the weed; s
m
a
x denotes the maximum number of seeds produced by each strain of weed; s
m
i
n denotes the minimum number of seeds produced per strain of weed; f
l
o
o
r
(
x
) 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 asN
s for each solution x
i in the parent generation, which will be added with normal random parameters, r
(
r
∈
[
-
σ
,
σ
]
), whose mean value μ
=
0 and the standard deviation is σ, respectively, so that each strain of weed will generate N
s 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 variabler is amenable to such a probability distribution whose mathematical expectation is μ and variance σ
2, denoting as N
(
μ
,
σ
2
), its probability density function can be described as follows:(21)
f
x
=
1
σ
2
π
e
x
p
-
x
-
μ
2
2
σ
2
.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 r 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 r
=
r
∗
σ
+
μ.The method of seed production is shown as follows:(22)
S
i
,
j
,
k
g
=
W
i
,
k
g
+
R
i
,
j
,
k
g
,where, i
=
1,2
,
…
,
P, g means the current generation. If g
=
1, P
=
P
m
i
n; otherwise P
=
P
m
a
x, j
=
1,2
,
…
,
N
s, k
=
1,2
,
…
,
d. W denotes weed space, S denotes seed space, and R 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:(23)
σ
=
G
-
g
G
n
σ
i
n
i
t
-
σ
f
i
n
a
l
+
σ
f
i
a
n
l
,where σ
i
n
i
t denotes the initial value of the interval step, σ
f
i
n
a
l denotes the final value of the interval step, σ denotes the current value of stepsize, G denotes the maximum generation time, g denotes the current generation, and n 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 populationP
m
a
x 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 P
m
a
x 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 f
(
x
) denote an objective function, and x
i denotes the group of points in the current simplex, i
∈
{
1
,
…
,
N
+
1
}. The NMSM method is inclusive of the following steps [34].Step 1 (order).
The system of inequalities,f
(
x
1
)
≤
f
(
x
2
)
≤
⋯
≤
f
(
x
i
)
≤
⋯
≤
f
(
x
N
+
1
), must meet the demands of the N
+
1 vertices.Step 2 (reflection).
Firstly, compute the reflection pointx
R of the simplex by x
R
=
x
¯
+
α
(
x
¯
-
x
N
+
1
), (
α
>
0
), where x
¯
=
(
1
/
N
)
∑
i
=
1
N is the centroid of the n best points. Then calculate f
R
=
f
(
x
R
). If f
1
≤
f
R
≤
f
N, accept the reflected point x
R and terminate the iteration.Step 3 (expandsion).
Iff
R
<
f
1, calculate the expansion point x
E of the simplex by x
E
=
x
¯
+
β
(
x
R
-
x
¯
), (β
>
1 and β
>
α). Compute f
E
=
f
(
x
E
). If f
E
<
f
R, accept x
E and terminate the iteration; otherwise, accept x
R and terminate the iteration.Step 4 (contraction).
Iff
R
≥
f
N, implement a contraction operation between x
¯ and the better point of x
R and x
n
+
1.(a)
Outside contraction: iff
N
≤
f
R
<
f
n
+
1, calculate the outside contraction point x
C by x
C
=
x
¯
+
γ
(
x
R
-
x
¯
), (0
<
γ
<
1). Compute f
C
=
f
(
x
C
). If f
C
<
f
R, accept x
C and terminate the iteration; otherwise, continue with Step 5.
(b)
Inside contraction: iff
R
≥
f
N
+
1, compute the inside contraction point x
C
C by x
C
C
=
x
¯
-
γ
(
x
¯
-
x
N
+
1
). Calculate f
C
C
=
f
(
x
C
C
). If f
C
C
<
f
R, accept x
C
C and end the iteration; if not, continue with Step 5.Step 5 (shrink).
Calculate theN points by v
i
=
x
1
+
σ
(
x
i
-
x
1
)
(
0
<
σ
<
1
). And calculate f
(
v
i
), i
=
2
,
…
,
N
+
1. The unordered vertices of the simplex at the next iteration consist of x
1
,
v
2
,
…
,
v
N
+
1. 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 Algorithm2, where W is the weed populations and b is the best individual.Algorithm 2 (HIWO algorithm).
RequireW, G, P
m
a
x, P
m
i
n, d, n, σ
i
n
i
t, σ
f
i
n
a
l, s
m
a
x, s
m
i
n, x
m
i
n, and x
m
a
x
.
Ensureb
.
(1)
Extract the set of reliable negative and/or positive samplesT
n from U
n with the help of P
n.
(2)
Initialize the weed populationW
0 in accordance with formula (19), let g
=
0.
(3)
Calculate the adaptability of weeds ofW
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).
(6)
Calculate the seed number produced by each individual of parent weedW
g using formula (20); compute the fitness of seeds S
g according to formula (22) until the numbers of W
g and S
g are greater than or equal to the maximum population P
m
a
x.
(7)
Compute the adaptability of seedsS
g.
(8)
After sortingW
g and S
g in accordance with the fitness value, choose the number P
m
a
x of the best results as the next generation weed W
g
+
1.
(9)
Update the best individualb
g of current generation.
(10)
Execute NMSM algorithm and implement a precise search.
(11)
Letg
=
g
+
1; 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]. F
1
–
F
9 belong to multipeak functions, and each of them owns multiple local extremum, in which F
7 has 100
!
=
9.33
E
+
157 local extreme points, and F
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 F
1 to F
20 in literature [35].Table 1
The characteristics of benchmark functions.
F
Domain
Optimum
D
F
Domain
Optimum
D
F
1
[
-
500,500
]
−12569.5
30
F
11
[
-
100,100
]
0
30
F
2
[
-
5.12,5.12
]
0
30
F
12
[
-
1.28,1.28
]
0
30
F
3
[
-
32,32
]
0
30
F
13
[
-
10,10
]
0
30
F
4
[
-
600,600
]
0
30
F
14
[
-
100,100
]
0
30
F
5
[
-
50,50
]
0
30
F
15
[
-
100,100
]
0
30
F
6
[
-
50,50
]
0
30
F
16
[
-
5,5
]
−1.031629
2
F
7
[
0
,
π
]
−99.2784
100
F
17
[
-
5,10
]
×
[
0,15
]
0.398
2
F
8
[
-
π
,
π
]
0
100
F
18
[
-
2,2
]
3
2
F
9
[
-
5,5
]
−78.33236
100
F
19
[
-
5,5
]
0.000308
4
F
10
[
-
5,10
]
0
100
F
20
[
0,1
]
−3.32
6The 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 F
1
–
F
15 are 1
E
5; the FES of IWO, PSO, and HIWO algorithms for functions F
16
–
F
20 are 1
E
4. The FES of the LEA algorithm for functions F
1
–
F
15 and F
16
–
F
20 are much more than 1
E
5 and 1
E
4, respectively.The parameter settings of the three intelligence algorithms (PSO, IWO, and HIWO) are as follows, whereP is the size of population, G is the maximum number of generations, and M
a
x
X is the maximum value of objective space. For the three algorithms, the G is the same, G
=
5000 (from F
1 to F
15) or G
=
500 (from F
16 to F
20); the setting of FES is the same.(i)
The PSO algorithm: population sizeP
=
20; learning factor C
1
=
C
2
=
1.49445; the initial weight value W
s
t
a
r
t
=
0.9; the final weight value W
e
n
d
=
0.4; the maximum velocity V
m
a
x
=
M
a
x
X
/
5.0; the minimum velocity V
m
i
n
=
-
M
a
x
X
/
5.0; the calculation formulation of weight factor for the current generation is W
=
W
s
t
a
r
t
-
(
W
s
t
a
r
t
-
W
e
n
d
)
∗
p
o
w
(
d
o
u
b
l
e
(
g
)
/
G
,
2
).
(ii)
The IWO algorithm: the minimum number of weed populationsP
m
i
n
=
8; the maximum number of weed populations P
m
a
x
=
10; the initial step Q
i
n
i
t
=
10; the final step Q
f
i
n
=
1
e
-
3; the minimum number of seeds S
m
i
n
=
2; the maximum number of seeds S
m
a
x
=
3; nonlinear index n
=
3.
(iii)
The HIWO algorithm: the minimum number of weed populationsP
m
i
n
=
4; the maximum number of weed populations P
m
a
x
=
5; the settings of other parameters of HIWO are the same as IWO.Table2 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 F
16
–
F
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 F
16 and F
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 F
3 and F
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 F
3, F
5, F
6, F
8, F
9, F
12, and F
19.Table 2
Mean value and standard deviation comparison of different algorithms for benchmark functions.
F
BFGS
NMSM
IWO
Mean
std
Mean
std
Mean
std
F
1
-
2198.02
2.26
E
+
03
-
6521.42
3.80
E
+
02
-
7340.4
6.33
E
+
02
F
2
3.13
E
+
02
3.75
E
+
01
4.39
E
+
01
8.02
5.82
E
+
01
1.53
E
+
01
F
3
1.98
E
+
01
5.08
E
-
01
1.01
E
+
01
1.03
3.80
E
-
03
2.37
E
-
04
F
4
4.61
E
+
02
5.52
E
+
01
3.34
E
+
01
1.23
E
+
01
8.62
E
-
03
7.21
E
-
03
F
5
3.52
E
+
08
9.55
E
+
07
8.93
5.39
3.46
E
-
03
1.86
E
-
02
F
6
6.27
E
+
08
2.00
E
+
08
2.58
E
+
03
3.70
E
+
03
3.70
E
-
04
1.97
E
-
03
F
7
-
31.84
1.86
E
+
01
-
50.03
4.38
-
45.91
2.83
F
8
4.15
E
+
07
1.50
E
+
07
1.03
E
+
07
2.34
E
+
06
9.13
E
+
05
3.53
E
+
05
F
9
-
33.97
2.71
-
47.24
9.07
-
67.53
1.46
F
10
6.43
E
+
06
7.64
E
+
05
1.50
E
+
05
4.35
E
+
04
1.40
E
+
02
6.14
E
+
01
F
11
4.33
E
+
04
1.32
E
+
04
2.36
E
+
03
7.10
E
+
02
2.68
E
-
05
2.94
E
-
06
F
12
4.61
1.58
1.64
E
-
01
8.19
E
-
02
1.55
E
-
02
6.64
E
-
03
F
13
1.74
E
+
11
3.49
E
+
11
9.06
2.8
2.29
E
-
02
1.53
E
-
03
F
14
1.55
E
+
05
1.14
E
+
05
2.05
E
+
04
1.35
E
+
04
4.02
E
-
04
1.77
E
-
04
F
15
8.32
E
+
01
1.09
E
+
01
1.84
E
+
01
2.64
2.28
E
-
03
1.80
E
-
04
F
16
27.55
4.24
E
+
01
-
1.03163
7.65
E
-
18
-
1.03163
2.49
E
-
09
F
17
8.96
7.04
0.397888
2.01
E
-
06
0.397913
3.85
E
-
12
F
18
1246.7
1.82
E
+
03
3
2.82
E
-
17
3
9.89
E
-
08
F
19
1.09
E
-
01
6.03
E
-
02
6.31
E
-
04
4.99
E
-
04
1.12
E
-
03
2.16
E
-
04
F
20
-
1.1488
4.12
E
-
01
-
3.2375
9.02
E
-
02
-
3.2031
2.15
E
-
06
F
PSO
LEA
HIWO
Mean
std
Mean
std
Mean
std
F
1
-
7307.54
6.00
E
+
02
-
12569.45
4.83
E
-
04
-
12569.49
1.68
E
-
13
F
2
4.29
E
+
01
1.20
E
+
01
2.10
E
-
18
3.36
E
-
18
0
0
F
3
1.14
E
-
01
3.47
E
-
01
3.27
E
-
16
3.00
E
-
17
2.37
E
-
15
1.67
E
-
16
F
4
1.43
E
-
02
1.68
E
-
02
6.10
E
-
16
2.51
E
-
17
0
0
F
5
6.91
E
-
03
2.59
E
-
02
2.48
E
-
06
2.28
E
-
06
8.03
E
-
17
1.16
E
-
27
F
6
5.49
E
-
03
1.13
E
-
02
1.73
E
-
04
1.21
E
-
04
6.33
E
-
19
4.48
E
-
27
F
7
-
69.32
4.61
-
93.01
2.31
E
-
02
-
90.38
1.86
E
-
02
F
8
6.23
E
+
06
3.41
E
+
06
1.63
E
-
06
6.53
E
-
07
3.28
E
-
13
3.56
E
-
13
F
9
-
66.07
1.34
-
78.31
6.13
E
-
03
-
78.3323
8.01
E
-
18
F
10
5.07
E
+
04
5.04
E
+
04
5.61
E
-
01
1.08
E
-
01
1.87
E
-
35
2.64
E
-
35
F
11
7.70
E
-
12
3.97
E
-
11
4.73
E
-
16
6.22
E
-
17
0
0
F
12
2.23
E
-
03
7.70
E
-
04
5.14
E
-
03
4.43
E
-
04
3.29
E
-
04
2.36
E
-
04
F
13
2.38
E
-
05
5.97
E
-
05
4.25
E
-
19
4.24
E
-
19
0
0
F
14
1.68
E
+
02
8.98
E
+
02
6.78
E
-
18
5.43
E
-
18
2.23
E
-
283
3.16
E
-
283
F
15
5.29
E
-
01
2.22
E
-
01
2.68
E
-
16
6.26
E
-
17
1.50
E
-
160
2.13
E
-
160
F
16
-
1.03163
3.42
E
-
18
-
1.03108
3.36
E
-
07
-
1.031628
1.01
E
-
18
F
17
0.397889
2.68
E
-
06
0.398
2.65
E
-
08
0.397887
0
F
18
3
1.05
E
-
17
3.00003
6.25
E
-
05
3
2.51
E
-
18
F
19
3.27
E
-
03
6.72
E
-
03
3.51
E
-
04
7.36
E
-
05
3.12
E
-
04
5.93
E
-
06
F
20
-
3.2451
8.59
E
-
02
-
3.301
7.83
E
-
03
-
3.322
1.38
E
-
17Table3 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.Table 3
Wilcoxon statistic result comparison by HIWO and other algorithms.
Metric
BFGS
NMSM
IWO
PSO
LEA
B
20
19
18
18
18
E
—
1
2
2
—
W
—
—
—
—
2It can be seen from Figure1 to Figure 20 that, compared with IWO and PSO, HIWO has faster convergence speed in solving benchmark functions.Figure 1
Convergence ofF
1.Figure 2
Convergence ofF
2.Figure 3
Convergence ofF
3.Figure 4
Convergence ofF
4.Figure 5
Convergence ofF
5.Figure 6
Convergence ofF
6.Figure 7
Convergence ofF
7.Figure 8
Convergence ofF
8.Figure 9
Convergence ofF
9.Figure 10
Convergence ofF
10.Figure 11
Convergence ofF
11.Figure 12
Convergence ofF
12.Figure 13
Convergence ofF
13.Figure 14
Convergence ofF
14.Figure 15
Convergence ofF
15.Figure 16
Convergence ofF
16.Figure 17
Convergence ofF
17.Figure 18
Convergence ofF
18.Figure 19
Convergence ofF
19.Figure 20
Convergence ofF
20.The convergence speed of HIWO is faster than those of IWO and PSO except forF
16 and F
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.


### 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Δ
t
=
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 θ. 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:(24)
A
A
E
=
1
n
∑
i
=
1
n
Q
~
t
i
-
Q
t
i
,
A
R
E
=
1
n
∑
i
=
1
n
Q
~
t
i
-
Q
t
i
Q
t
i
.It can be seen from Table 4 that the experiment results are much better, when θ
=
1
/
3, confirming the error estimation in this paper.Table 4
Numerical results with differentθ for flood routing in 1960.
θ
AAE
ARE (%)
θ
AAE
ARE (%)
θ
=
0
6.844292
1.964226
θ
=
1
/
7
6.813278
2.010455
θ
=
1
/
6
19.629630
5.454065
θ
=
2
/
3
6.838390
2.007756
θ
=
1
/
5
7.037453
2.039629
θ
=
3
/
4
6.837713
2.006820
θ
=
1
/
4
7.108511
2.090138
θ
=
5
/
6
6.843642
2.006165
θ
=
1
/
3
6.829915
2.010007
θ
=
1
6.833480
2.004863On 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 [
K
,
x
,
m
]
=
[
0,1
,
1
] and θ
=
1
/
3. Once the parameters K, x, and m are calculated, we can use (16) to calculate the next calculated outflow Q
~
(
t
i
), where Q
~
(
t
1
)
=
Q
(
t
1
). Obviously, given that (16) is a nonlinear equation about Q
(
t
i
), 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.Table 5
Result comparison of different methods in 1960.
Algorithm
K
x
m
AAE
ARE (%)
BFGS [5]
10.815576
0.982641
1.021860
7.866438
2.345940
PSO [6]
13.095689
0.800769
0.995361
6.831708
2.010701
NMSM [9]
0.001628
0.758777
2.337256
9.871414
3.308699
IWO
4.998601
0.881355
1.136204
7.241756
2.148295
HIWO
13.281324
0.800085
0.993263
6.829915
2.010007Table 6
Result comparison of different methods in 1961 and 1964.
Algorithm
1961
1964
AAE
ARE (%)
AAE
ARE (%)
BFGS [5]
4.634922
1.156255
8.745983
2.185487
PSO [6]
3.980577
0.991698
8.710431
2.179195
NMSM [9]
10.990222
2.808921
5.334347
1.709903
IWO
4.154933
1.064907
8.475477
2.136220
HIWO
3.971981
0.990858
8.716951
2.180527As the numerical results in Tables5 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 Figures21, 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.Figure 21
Comparison between calculated values and observed values in 1960.Figure 22
Comparison between calculated values and observed values in 1961.Figure 23
Comparison between calculated values and observed values in 1964.


## 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]. F
1
–
F
9 belong to multipeak functions, and each of them owns multiple local extremum, in which F
7 has 100
!
=
9.33
E
+
157 local extreme points, and F
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 F
1 to F
20 in literature [35].Table 1
The characteristics of benchmark functions.
F
Domain
Optimum
D
F
Domain
Optimum
D
F
1
[
-
500,500
]
−12569.5
30
F
11
[
-
100,100
]
0
30
F
2
[
-
5.12,5.12
]
0
30
F
12
[
-
1.28,1.28
]
0
30
F
3
[
-
32,32
]
0
30
F
13
[
-
10,10
]
0
30
F
4
[
-
600,600
]
0
30
F
14
[
-
100,100
]
0
30
F
5
[
-
50,50
]
0
30
F
15
[
-
100,100
]
0
30
F
6
[
-
50,50
]
0
30
F
16
[
-
5,5
]
−1.031629
2
F
7
[
0
,
π
]
−99.2784
100
F
17
[
-
5,10
]
×
[
0,15
]
0.398
2
F
8
[
-
π
,
π
]
0
100
F
18
[
-
2,2
]
3
2
F
9
[
-
5,5
]
−78.33236
100
F
19
[
-
5,5
]
0.000308
4
F
10
[
-
5,10
]
0
100
F
20
[
0,1
]
−3.32
6The 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 F
1
–
F
15 are 1
E
5; the FES of IWO, PSO, and HIWO algorithms for functions F
16
–
F
20 are 1
E
4. The FES of the LEA algorithm for functions F
1
–
F
15 and F
16
–
F
20 are much more than 1
E
5 and 1
E
4, respectively.The parameter settings of the three intelligence algorithms (PSO, IWO, and HIWO) are as follows, whereP is the size of population, G is the maximum number of generations, and M
a
x
X is the maximum value of objective space. For the three algorithms, the G is the same, G
=
5000 (from F
1 to F
15) or G
=
500 (from F
16 to F
20); the setting of FES is the same.(i)
The PSO algorithm: population sizeP
=
20; learning factor C
1
=
C
2
=
1.49445; the initial weight value W
s
t
a
r
t
=
0.9; the final weight value W
e
n
d
=
0.4; the maximum velocity V
m
a
x
=
M
a
x
X
/
5.0; the minimum velocity V
m
i
n
=
-
M
a
x
X
/
5.0; the calculation formulation of weight factor for the current generation is W
=
W
s
t
a
r
t
-
(
W
s
t
a
r
t
-
W
e
n
d
)
∗
p
o
w
(
d
o
u
b
l
e
(
g
)
/
G
,
2
).
(ii)
The IWO algorithm: the minimum number of weed populationsP
m
i
n
=
8; the maximum number of weed populations P
m
a
x
=
10; the initial step Q
i
n
i
t
=
10; the final step Q
f
i
n
=
1
e
-
3; the minimum number of seeds S
m
i
n
=
2; the maximum number of seeds S
m
a
x
=
3; nonlinear index n
=
3.
(iii)
The HIWO algorithm: the minimum number of weed populationsP
m
i
n
=
4; the maximum number of weed populations P
m
a
x
=
5; the settings of other parameters of HIWO are the same as IWO.Table2 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 F
16
–
F
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 F
16 and F
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 F
3 and F
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 F
3, F
5, F
6, F
8, F
9, F
12, and F
19.Table 2
Mean value and standard deviation comparison of different algorithms for benchmark functions.
F
BFGS
NMSM
IWO
Mean
std
Mean
std
Mean
std
F
1
-
2198.02
2.26
E
+
03
-
6521.42
3.80
E
+
02
-
7340.4
6.33
E
+
02
F
2
3.13
E
+
02
3.75
E
+
01
4.39
E
+
01
8.02
5.82
E
+
01
1.53
E
+
01
F
3
1.98
E
+
01
5.08
E
-
01
1.01
E
+
01
1.03
3.80
E
-
03
2.37
E
-
04
F
4
4.61
E
+
02
5.52
E
+
01
3.34
E
+
01
1.23
E
+
01
8.62
E
-
03
7.21
E
-
03
F
5
3.52
E
+
08
9.55
E
+
07
8.93
5.39
3.46
E
-
03
1.86
E
-
02
F
6
6.27
E
+
08
2.00
E
+
08
2.58
E
+
03
3.70
E
+
03
3.70
E
-
04
1.97
E
-
03
F
7
-
31.84
1.86
E
+
01
-
50.03
4.38
-
45.91
2.83
F
8
4.15
E
+
07
1.50
E
+
07
1.03
E
+
07
2.34
E
+
06
9.13
E
+
05
3.53
E
+
05
F
9
-
33.97
2.71
-
47.24
9.07
-
67.53
1.46
F
10
6.43
E
+
06
7.64
E
+
05
1.50
E
+
05
4.35
E
+
04
1.40
E
+
02
6.14
E
+
01
F
11
4.33
E
+
04
1.32
E
+
04
2.36
E
+
03
7.10
E
+
02
2.68
E
-
05
2.94
E
-
06
F
12
4.61
1.58
1.64
E
-
01
8.19
E
-
02
1.55
E
-
02
6.64
E
-
03
F
13
1.74
E
+
11
3.49
E
+
11
9.06
2.8
2.29
E
-
02
1.53
E
-
03
F
14
1.55
E
+
05
1.14
E
+
05
2.05
E
+
04
1.35
E
+
04
4.02
E
-
04
1.77
E
-
04
F
15
8.32
E
+
01
1.09
E
+
01
1.84
E
+
01
2.64
2.28
E
-
03
1.80
E
-
04
F
16
27.55
4.24
E
+
01
-
1.03163
7.65
E
-
18
-
1.03163
2.49
E
-
09
F
17
8.96
7.04
0.397888
2.01
E
-
06
0.397913
3.85
E
-
12
F
18
1246.7
1.82
E
+
03
3
2.82
E
-
17
3
9.89
E
-
08
F
19
1.09
E
-
01
6.03
E
-
02
6.31
E
-
04
4.99
E
-
04
1.12
E
-
03
2.16
E
-
04
F
20
-
1.1488
4.12
E
-
01
-
3.2375
9.02
E
-
02
-
3.2031
2.15
E
-
06
F
PSO
LEA
HIWO
Mean
std
Mean
std
Mean
std
F
1
-
7307.54
6.00
E
+
02
-
12569.45
4.83
E
-
04
-
12569.49
1.68
E
-
13
F
2
4.29
E
+
01
1.20
E
+
01
2.10
E
-
18
3.36
E
-
18
0
0
F
3
1.14
E
-
01
3.47
E
-
01
3.27
E
-
16
3.00
E
-
17
2.37
E
-
15
1.67
E
-
16
F
4
1.43
E
-
02
1.68
E
-
02
6.10
E
-
16
2.51
E
-
17
0
0
F
5
6.91
E
-
03
2.59
E
-
02
2.48
E
-
06
2.28
E
-
06
8.03
E
-
17
1.16
E
-
27
F
6
5.49
E
-
03
1.13
E
-
02
1.73
E
-
04
1.21
E
-
04
6.33
E
-
19
4.48
E
-
27
F
7
-
69.32
4.61
-
93.01
2.31
E
-
02
-
90.38
1.86
E
-
02
F
8
6.23
E
+
06
3.41
E
+
06
1.63
E
-
06
6.53
E
-
07
3.28
E
-
13
3.56
E
-
13
F
9
-
66.07
1.34
-
78.31
6.13
E
-
03
-
78.3323
8.01
E
-
18
F
10
5.07
E
+
04
5.04
E
+
04
5.61
E
-
01
1.08
E
-
01
1.87
E
-
35
2.64
E
-
35
F
11
7.70
E
-
12
3.97
E
-
11
4.73
E
-
16
6.22
E
-
17
0
0
F
12
2.23
E
-
03
7.70
E
-
04
5.14
E
-
03
4.43
E
-
04
3.29
E
-
04
2.36
E
-
04
F
13
2.38
E
-
05
5.97
E
-
05
4.25
E
-
19
4.24
E
-
19
0
0
F
14
1.68
E
+
02
8.98
E
+
02
6.78
E
-
18
5.43
E
-
18
2.23
E
-
283
3.16
E
-
283
F
15
5.29
E
-
01
2.22
E
-
01
2.68
E
-
16
6.26
E
-
17
1.50
E
-
160
2.13
E
-
160
F
16
-
1.03163
3.42
E
-
18
-
1.03108
3.36
E
-
07
-
1.031628
1.01
E
-
18
F
17
0.397889
2.68
E
-
06
0.398
2.65
E
-
08
0.397887
0
F
18
3
1.05
E
-
17
3.00003
6.25
E
-
05
3
2.51
E
-
18
F
19
3.27
E
-
03
6.72
E
-
03
3.51
E
-
04
7.36
E
-
05
3.12
E
-
04
5.93
E
-
06
F
20
-
3.2451
8.59
E
-
02
-
3.301
7.83
E
-
03
-
3.322
1.38
E
-
17Table3 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.Table 3
Wilcoxon statistic result comparison by HIWO and other algorithms.
Metric
BFGS
NMSM
IWO
PSO
LEA
B
20
19
18
18
18
E
—
1
2
2
—
W
—
—
—
—
2It can be seen from Figure1 to Figure 20 that, compared with IWO and PSO, HIWO has faster convergence speed in solving benchmark functions.Figure 1
Convergence ofF
1.Figure 2
Convergence ofF
2.Figure 3
Convergence ofF
3.Figure 4
Convergence ofF
4.Figure 5
Convergence ofF
5.Figure 6
Convergence ofF
6.Figure 7
Convergence ofF
7.Figure 8
Convergence ofF
8.Figure 9
Convergence ofF
9.Figure 10
Convergence ofF
10.Figure 11
Convergence ofF
11.Figure 12
Convergence ofF
12.Figure 13
Convergence ofF
13.Figure 14
Convergence ofF
14.Figure 15
Convergence ofF
15.Figure 16
Convergence ofF
16.Figure 17
Convergence ofF
17.Figure 18
Convergence ofF
18.Figure 19
Convergence ofF
19.Figure 20
Convergence ofF
20.The convergence speed of HIWO is faster than those of IWO and PSO except forF
16 and F
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.


## 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Δ
t
=
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 θ. 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:(24)
A
A
E
=
1
n
∑
i
=
1
n
Q
~
t
i
-
Q
t
i
,
A
R
E
=
1
n
∑
i
=
1
n
Q
~
t
i
-
Q
t
i
Q
t
i
.It can be seen from Table 4 that the experiment results are much better, when θ
=
1
/
3, confirming the error estimation in this paper.Table 4
Numerical results with differentθ for flood routing in 1960.
θ
AAE
ARE (%)
θ
AAE
ARE (%)
θ
=
0
6.844292
1.964226
θ
=
1
/
7
6.813278
2.010455
θ
=
1
/
6
19.629630
5.454065
θ
=
2
/
3
6.838390
2.007756
θ
=
1
/
5
7.037453
2.039629
θ
=
3
/
4
6.837713
2.006820
θ
=
1
/
4
7.108511
2.090138
θ
=
5
/
6
6.843642
2.006165
θ
=
1
/
3
6.829915
2.010007
θ
=
1
6.833480
2.004863On 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 [
K
,
x
,
m
]
=
[
0,1
,
1
] and θ
=
1
/
3. Once the parameters K, x, and m are calculated, we can use (16) to calculate the next calculated outflow Q
~
(
t
i
), where Q
~
(
t
1
)
=
Q
(
t
1
). Obviously, given that (16) is a nonlinear equation about Q
(
t
i
), 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.Table 5
Result comparison of different methods in 1960.
Algorithm
K
x
m
AAE
ARE (%)
BFGS [5]
10.815576
0.982641
1.021860
7.866438
2.345940
PSO [6]
13.095689
0.800769
0.995361
6.831708
2.010701
NMSM [9]
0.001628
0.758777
2.337256
9.871414
3.308699
IWO
4.998601
0.881355
1.136204
7.241756
2.148295
HIWO
13.281324
0.800085
0.993263
6.829915
2.010007Table 6
Result comparison of different methods in 1961 and 1964.
Algorithm
1961
1964
AAE
ARE (%)
AAE
ARE (%)
BFGS [5]
4.634922
1.156255
8.745983
2.185487
PSO [6]
3.980577
0.991698
8.710431
2.179195
NMSM [9]
10.990222
2.808921
5.334347
1.709903
IWO
4.154933
1.064907
8.475477
2.136220
HIWO
3.971981
0.990858
8.716951
2.180527As the numerical results in Tables5 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 Figures21, 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.Figure 21
Comparison between calculated values and observed values in 1960.Figure 22
Comparison between calculated values and observed values in 1961.Figure 23
Comparison between calculated values and observed values in 1964.


## 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 parametersK, x, and m 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.



---

*Source: 573894-2015-07-02.xml*