# A Novel Parameter Estimation Method for Muskingum Model Using New Newton-Type Trust Region Algorithm

Authors: Zhou Sheng; Aijia Ouyang; Li-Bin Liu; Gonglin Yuan
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Parameters estimation of Muskingum model is very significative in both exploitation and utilization of water resources and hydrological forecasting. The optimal results of parameters directly affect the accuracy of flood forecasting. This paper considers the parameters estimation problem of Muskingum model from the following two aspects. Firstly, based on the general trapezoid formulas, a class of new discretization methods including a parameterθ to approximate Muskingum model is presented. The accuracy of these methods is second-order, when θ
≠
1
/
3. Particularly, if we choose θ
=
1
/
3, the accuracy of the presented method can be improved to third-order. Secondly, according to the Newton-type trust region algorithm, a new Newton-type trust region algorithm is given to obtain the parameters of Muskingum model. This method can avoid high dependence on the initial parameters. The average absolute errors (AAE) and the average relative errors (ARE) of the proposed algorithm of parameters estimation for Muskingum model are 8.208122 and 2.462438%, respectively, where θ
=
1
/
3. It is shown from these results that the presented algorithm has higher forecasting accuracy and wider practicability than other methods.


---

## Body

## 1. Introduction

Flood routing in open channels is a very important tool in the design of flood protection measures to estimate how the proposed measures will affect the behavior of flood waves in rivers so that enough protection and economic solutions can be found [1]. In general, flood routing procedures can be classified as either hydrologic or hydraulic. Hydrologic routing method is on the basis of the storage-continuity equation, whereas hydraulic routing method is on the basis of both continuity and momentum equations. One of the hydrologic routing approaches, the Muskingum method, was first developed by McCarthy for flood control studies in the Muskingum river basin in Ohio.The following continuity and storage equations are used to describe the Muskingum model:(1)
d
W
(
t
)
d
t
=
I
(
t
)
-
Q
(
t
)
,
where W
(
t
) represents the channel storage at time t; I
(
t
) and Q
(
t
) represent the rates of inflow and outflow at time t, respectively. The linear Musikingum model is
(2)
W
(
t
)
=
K
x
I
t
+
1
-
x
Q
t
,
where K is a storage time constant for the river reach and x is a weighting factor commonly varying between 0.0 and 0.3 for the river channel.It is worth mentioning that the parametersK and x in (2) are graphically estimated by a trial-and-error procedure [2]. 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 and plotted against W. The particular value which generates the loop is accepted as the best estimation of x. The slope of the straight line fits through the loop derives K. Although the trial-and-error procedure has been used for many years, it is time consuming and subjective interpretation. Therefore, in order to avoid subjective interpretations of observed data in estimating K and x, the following finite difference scheme [3] is used to the resulting ordinary differential equation (1):
(3)
W
(
t
i
)
-
W
(
t
i
-
1
)
Δ
t
=
I
(
t
i
)
+
i
(
t
i
-
1
)
2
-
Q
(
t
i
)
+
Q
(
t
i
-
1
)
2
,
where
(4)
W
(
t
i
)
=
K
[
x
I
(
t
i
)
+
(
1
-
x
)
Q
(
t
i
)
]
,
W
t
i
-
1
=
K
x
I
t
i
-
1
+
1
-
x
Q
t
i
-
1
,
i
=
2
,
…
,
n
,
and Q
(
t
i
) represents observed outflow discharges at time t
i, I
(
t
i
) represents the inflow discharge at time interval t
i, Δ
t is the time step, and n is the total time number. Substituting (4) into (3), we have
(5)
Q
(
t
i
+
1
)
=
c
0
I
(
t
i
+
1
)
+
c
1
I
(
t
i
)
+
c
2
Q
(
t
i
)
,
where
(6)
c
0
=
-
K
x
+
Δ
t
/
2
K
(
1
-
x
)
+
Δ
t
/
2
,
c
1
=
K
x
+
Δ
t
/
2
K
(
1
-
x
)
+
Δ
t
/
2
,
c
2
=
K
(
1
-
x
)
-
Δ
t
/
2
K
(
1
-
x
)
+
Δ
t
/
2
.By minimizing the sum of the square of the deviations between observed and calculated outflows, we obtain the following objective function:(7)
min
⁡
f
=
∑
i
=
1
n
-
1
‍
Q
(
t
i
+
1
)
-
c
0
I
(
t
i
+
1
)
-
c
1
I
(
t
i
)
-
c
2
Q
(
t
i
)
2
.In the past two decades, in order to obtain the parametersK and x, some researchers adopted many optimization methods to solve the above optimization problem (7). Until now, these methods can be classified into traditional optimization methods and intelligent algorithms. The existing traditional methods include nonlinear programming method (NPM) [4], the least-square method (L-SM) [5, 6], method of trial-and-error (TAE) [2], the minimum area method (MAM) [1, 7], the Broydene-Fletchere-Goldfarbe-Shanno (BFGS) technique [8], test-method and least residual square method [9], and Nelder-Mead simplex method [10]. These traditional methods have their own advantages, but there also exist complex calculations and poor generality disadvantages, and some of these methods are related to the selection of initial point, which is easy to fall into local optimum. The existing intelligent algorithms include genetic algorithm [11], harmony search [12], Gray-encoded accelerating genetic algorithm [9], particle swarm optimization [13], immune clonal selection algorithm [14], and differential evolution algorithm [15]. The advantages of those intelligent algorithms are that the global search ability is strong and the program is relatively simple to design. However, these methods have unavoidable disadvantages of slow convergence precision, low precision solution in limited generations, and premature convergence.In a word, almost all of the above methods used to estimate parametersK and x are based on (3). We know that the truncation error of (3) is only O
(
Δ
t
)
2. Thus, it is very necessary to design a higher accuracy method to discrete differential equation (1). Based on the generalized trapezoidal formula [16], this paper will first develop a class of new difference scheme which contains a parameter θ to approximate the differential equation (1). Then, similar to the above optimization problem (7), we can obtain a class of new unconstrained nonlinear optimization problems which contain parameter θ. It is noted that the accuracy of the presented difference schemes to approximate the differential equation (1) can be improved to third order, when θ
=
1
/
3. In other words, we can get a higher accuracy parameter estimation model, if θ
=
1
/
3. In addition, based on the Newton-type trust region algorithm [17], this paper will design a new Newton-type trust region algorithm (NN-TTRA) for estimating the Muskingum model parameters. Briefly speaking, this algorithm can avoid high dependence on the initial parameters and find the global optimization solution quickly.


## 2. A New Parameter Estimation of Muskingum Model

### 2.1. Model Description

At first, (1) can be rewritten as follows:
(8)
d
W
(
t
)
d
t
=
W
(
t
)
+
P
(
t
)
,
where P
(
t
)
=
I
(
t
)
-
Q
(
t
)
-
W
(
t
).Then, applying the generalized trapezoidal formula [16] GTF (θ) of Chawla et al. to (8), we can obtain
(9)
W
~
(
t
i
)
=
W
(
t
i
+
1
)
-
Δ
t
d
W
(
t
i
+
1
)
d
t
,
(10)
W
t
i
+
1
=
W
(
t
i
)
+
Δ
t
2
1
-
θ
d
W
t
i
d
t
hhhhhhhhhhhh
+
θ
d
W
~
(
t
i
)
d
t
+
d
W
(
t
i
+
1
)
d
t
.
Combining (8), (9), and (10), we get
(11)
W
~
(
t
i
)
=
(
1
-
Δ
t
)
W
(
t
i
+
1
)
-
Δ
t
P
(
t
i
+
1
)
,
(12)
W
t
i
+
1
=
W
(
t
i
)
+
Δ
t
2
P
(
t
i
+
1
)
+
P
(
t
i
)
+
(
1
-
θ
)
W
(
t
i
)
W
~
h
h
h
h
h
h
h
h
h
h
h
h
+
θ
W
~
(
t
i
)
+
W
(
t
i
+
1
)
.
Substituting (11) into (12), we finally obtain
(13)
1
-
Δ
t
2
θ
W
t
i
+
1
=
1
-
Δ
t
2
θ
W
(
t
i
)
+
Δ
t
2
I
(
t
i
)
-
Q
(
t
i
)
+
Δ
t
2
1
-
Δ
t
θ
I
(
t
i
+
1
)
-
Q
(
t
i
+
1
)
.
Obviously, when θ
=
0, (13) becomes (3).At last, substituting (4) into (13), we have
(14)
Q
(
t
i
+
1
)
=
a
0
I
(
t
i
+
1
)
+
a
1
I
(
t
i
)
+
a
2
Q
(
t
i
)
,
where
(15)
a
0
=
-
K
x
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
,
a
1
=
K
x
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
,
a
2
=
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
-
Δ
t
/
2
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
.By minimizing the residual sum of squares between observed and calculated outflows, the objective function can be given as follows:(16)
min
⁡
f
=
∑
i
=
1
n
-
1
‍
Q
(
t
i
+
1
)
-
a
0
I
(
t
i
+
1
)
-
a
1
I
(
t
i
)
-
a
2
Q
(
t
i
)
2
.


### 2.2. Truncation Errors Analysis

Here, we analyze the local truncation error of (13) at time direction [16]. Firstly, from (9), we have
(17)
d
W
~
(
t
i
)
d
t
=
d
W
t
i
+
1
d
t
-
Δ
t
d
2
W
t
i
+
1
d
t
2
.
Then, by using Taylor expansion, (10) can be written as follows:
(18)
-
Δ
t
3
12
d
3
W
t
i
d
t
4
-
Δ
t
4
24
d
4
W
t
i
d
t
4
+
O
Δ
t
5
=
Δ
t
θ
2
d
W
~
t
i
d
t
-
d
W
t
i
d
t
+
O
Δ
t
5
.
Substituting for d
W
~
(
t
i
)
/
d
t from (17) into (18), the truncation errors of generalized trapezoidal formula can be obtained as follows:
(19)
E
=
θ
4
-
1
12
Δ
t
3
d
3
W
(
t
i
)
d
t
3
+
θ
6
-
1
24
Δ
t
4
d
4
W
(
t
i
)
d
t
4
+
O
Δ
t
5
.From (19), it is clear that the order of our presented method in (13) is O
(
Δ
t
)
2 if θ
≠
1
/
3. In particular, for θ
=
1
/
3, our presented method in (13) is third-order accuracy.


## 2.1. Model Description

At first, (1) can be rewritten as follows:
(8)
d
W
(
t
)
d
t
=
W
(
t
)
+
P
(
t
)
,
where P
(
t
)
=
I
(
t
)
-
Q
(
t
)
-
W
(
t
).Then, applying the generalized trapezoidal formula [16] GTF (θ) of Chawla et al. to (8), we can obtain
(9)
W
~
(
t
i
)
=
W
(
t
i
+
1
)
-
Δ
t
d
W
(
t
i
+
1
)
d
t
,
(10)
W
t
i
+
1
=
W
(
t
i
)
+
Δ
t
2
1
-
θ
d
W
t
i
d
t
hhhhhhhhhhhh
+
θ
d
W
~
(
t
i
)
d
t
+
d
W
(
t
i
+
1
)
d
t
.
Combining (8), (9), and (10), we get
(11)
W
~
(
t
i
)
=
(
1
-
Δ
t
)
W
(
t
i
+
1
)
-
Δ
t
P
(
t
i
+
1
)
,
(12)
W
t
i
+
1
=
W
(
t
i
)
+
Δ
t
2
P
(
t
i
+
1
)
+
P
(
t
i
)
+
(
1
-
θ
)
W
(
t
i
)
W
~
h
h
h
h
h
h
h
h
h
h
h
h
+
θ
W
~
(
t
i
)
+
W
(
t
i
+
1
)
.
Substituting (11) into (12), we finally obtain
(13)
1
-
Δ
t
2
θ
W
t
i
+
1
=
1
-
Δ
t
2
θ
W
(
t
i
)
+
Δ
t
2
I
(
t
i
)
-
Q
(
t
i
)
+
Δ
t
2
1
-
Δ
t
θ
I
(
t
i
+
1
)
-
Q
(
t
i
+
1
)
.
Obviously, when θ
=
0, (13) becomes (3).At last, substituting (4) into (13), we have
(14)
Q
(
t
i
+
1
)
=
a
0
I
(
t
i
+
1
)
+
a
1
I
(
t
i
)
+
a
2
Q
(
t
i
)
,
where
(15)
a
0
=
-
K
x
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
,
a
1
=
K
x
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
,
a
2
=
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
-
Δ
t
/
2
K
(
1
-
x
)
(
1
-
Δ
t
/
2
θ
)
+
Δ
t
/
2
(
1
-
Δ
t
θ
)
.By minimizing the residual sum of squares between observed and calculated outflows, the objective function can be given as follows:(16)
min
⁡
f
=
∑
i
=
1
n
-
1
‍
Q
(
t
i
+
1
)
-
a
0
I
(
t
i
+
1
)
-
a
1
I
(
t
i
)
-
a
2
Q
(
t
i
)
2
.


## 2.2. Truncation Errors Analysis

Here, we analyze the local truncation error of (13) at time direction [16]. Firstly, from (9), we have
(17)
d
W
~
(
t
i
)
d
t
=
d
W
t
i
+
1
d
t
-
Δ
t
d
2
W
t
i
+
1
d
t
2
.
Then, by using Taylor expansion, (10) can be written as follows:
(18)
-
Δ
t
3
12
d
3
W
t
i
d
t
4
-
Δ
t
4
24
d
4
W
t
i
d
t
4
+
O
Δ
t
5
=
Δ
t
θ
2
d
W
~
t
i
d
t
-
d
W
t
i
d
t
+
O
Δ
t
5
.
Substituting for d
W
~
(
t
i
)
/
d
t from (17) into (18), the truncation errors of generalized trapezoidal formula can be obtained as follows:
(19)
E
=
θ
4
-
1
12
Δ
t
3
d
3
W
(
t
i
)
d
t
3
+
θ
6
-
1
24
Δ
t
4
d
4
W
(
t
i
)
d
t
4
+
O
Δ
t
5
.From (19), it is clear that the order of our presented method in (13) is O
(
Δ
t
)
2 if θ
≠
1
/
3. In particular, for θ
=
1
/
3, our presented method in (13) is third-order accuracy.


## 3. Algorithm Construction

### 3.1. Newton-Type Trust Region Algorithm

As far as we know, the basic trust region algorithm used to solve the following unconstrained optimization problem(20)
min
⁡
x
∈
R
n
f
(
x
)
was first presented clearly in [19]. Recently, many researchers have proposed some improved trust region algorithm to solve (20); see for example Esmaeili and Kimiaei [20] and Amini and Ahookhosh [21].Here, we assume thatx
j is the jth iterative point, f
j
=
f
(
x
j
), g
j
=
∇
f
(
x
j
), and B
j is the jth iteration of Hesse matrix ∇
2
f
(
x
j
); the trust region subproblem, which is at the jth iterative step of problem (20), can be formulated as follows:
(21)
min
⁡
q
j
(
d
)
=
g
j
T
d
+
1
2
d
T
B
j
d
,
.
s
.
t
.
.
.
d
≤
Δ
j
,
where Δ
j is a trust region radius and d
∈
R
n is an iterative step and · denotes the Euclidian norm of vectors or its induced matrix norm.Denote(22)
Δ
f
j
=
f
(
x
j
)
-
f
(
x
j
+
d
j
)
,
(23)
Δ
q
j
=
q
j
(
0
)
-
q
j
(
d
j
)
,
(24)
r
j
=
Δ
f
j
Δ
q
j
;
the general trust region algorithm [17] for solving the unconstrained optimization problem (20) is given as follows.Algorithm 1.
Step
0. Given a starting point x
0
∈
R
n, Δ
0
>
0 is the initial trust region radium, ɛ
>
0, B
0
=
Hesse
(
x
0
). Set j
:
=
0.
Step
1. Calculate g
j, if g
j
≤
ɛ, and stop iteration; otherwise, go to Step 2.
Step
2. Utilize the smoothing Newton method to obtain the solution d
j of the subproblem (21).
Step
3. Calculate the r
j of formula (24).
Step
4. Regulate trust region radius. If r
j
<
0.25, let Δ
j
+
1
:
=
0.5
Δ
j; if r
j
>
0.75 and d
j
=
Δ
j, let Δ
j
+
1
:
=
2
Δ
j; otherwise, let Δ
j
+
1
:
=
Δ
j.
Step
5. If r
j
>
0.25, let x
j
+
1
:
=
x
j
+
d
j, update B
j to be B
j
+
1
=
Hesse
(
x
j
+
1
), let j
:
=
j
+
1, and go to Step 1; otherwise, let x
j
+
1
=
x
j, j
:
=
j
+
1, and go to Step 2, where Δ
0
=
1
/
10
g
(
x
0
), ɛ
=
1
0
-
6.Obviously, the above trust region subproblem (21) is an unconstrained optimization problem whose objective function is a quadratic. For the sake of convenience, we first denote
(25)
z
=
μ
,
λ
,
d
,
(26)
Y
(
z
)
:
=
μ
λ
-
d
2
2
+
Δ
j
2
-
λ
+
d
2
2
-
Δ
j
2
2
+
4
μ
2
B
j
+
λ
E
d
-
g
j
,
(27)
β
(
z
)
=
γ
Y
z
min
⁡
1
,
Y
z
,
respectively. Then, the following smoothing Newton algorithm [17] is given to solve the subproblem (21).Algorithm 2.
Step
0. Given d
0
∈
R
n, z
0
=
(
μ
0
,
λ
0
,
d
0
), z
~
0
=
(
μ
0
,
0,0
). Select parameters δ
,
σ
,
γ
∈
(
0,1
), μ
0
γ
<
1, and γ
Y
(
z
0
)
<
1; set h
:
=
0.
Step
1. Calculate Y
(
z
h
), if Y
(
z
h
)
=
0, and terminate algorithm; otherwise, calculate β
h
=
β
(
z
h
).
Step
2. Obtain the solution for equations Y
(
z
h
)
+
Y
′
(
z
h
)
Δ
z
h
=
β
h
z
~, and the solution is Δ
z
h
=
(
Δ
μ
h
,
Δ
λ
h
,
Δ
d
h
).
Step
3. Suppose that m
h is the minimum nonnegative integer meeting Y
(
z
h
+
δ
m
h
Δ
z
h
)
≤
[
1
-
σ
(
1
-
β
μ
0
)
δ
m
h
]
Y
(
z
h
). Let α
h
:
=
δ
m
h and z
h
+
1
=
z
h
+
α
h
Δ
z
h.
Step
4. Let h
:
=
h
+
1, and go to Step 1.


### 3.2. A New Newton-Type Trust Region Algorithm

Optimizing ability of Newton-type trust region algorithm (N-TTRA) highly depends on initial parameters. And the algorithm is short of global optimizing capability although its local optimizing is fast. In order to enhance the global optimization capability and get rid of dependence on initial parameters, we construct a new Newton-type trust region algorithm by combining a new BFGS updating formula with N-TTRA. The basic idea is to updateB
j
+
1 with new BFGS formulas at Step 5 in Algorithm 1. This idea guarantees positive definiteness of B
j
+
1 to improve the global optimization capability; meanwhile, it gets rid of dependence on the initial parameter selection. New adjustment formula of BFGS is as follows:
(28)
B
j
+
1
=
B
j
+
y
j
*
y
j
*
T
s
j
T
y
j
*
-
B
j
s
j
s
j
T
B
j
s
j
T
B
j
s
j
,
where y
j
*
=
y
j
T
s
j
/
y
j
T
s
j
y
j, s
j
=
x
j
+
1
-
x
j, and y
j
=
g
j
+
1
-
g
j. Here, y
j
* guarantees positive definiteness of B
j
+
1, so the new Newton-type trust region algorithm (NN-TTRA) uses (28) as the new adjustment formula of B
j
+
1 at Step 5 in Algorithm 1. It is proved that the method is of super-linear convergence [22]. The flow chart of NN-TTRA is given in Figure 1.Figure 1
Flow chart of NN-TTRA.


## 3.1. Newton-Type Trust Region Algorithm

As far as we know, the basic trust region algorithm used to solve the following unconstrained optimization problem(20)
min
⁡
x
∈
R
n
f
(
x
)
was first presented clearly in [19]. Recently, many researchers have proposed some improved trust region algorithm to solve (20); see for example Esmaeili and Kimiaei [20] and Amini and Ahookhosh [21].Here, we assume thatx
j is the jth iterative point, f
j
=
f
(
x
j
), g
j
=
∇
f
(
x
j
), and B
j is the jth iteration of Hesse matrix ∇
2
f
(
x
j
); the trust region subproblem, which is at the jth iterative step of problem (20), can be formulated as follows:
(21)
min
⁡
q
j
(
d
)
=
g
j
T
d
+
1
2
d
T
B
j
d
,
.
s
.
t
.
.
.
d
≤
Δ
j
,
where Δ
j is a trust region radius and d
∈
R
n is an iterative step and · denotes the Euclidian norm of vectors or its induced matrix norm.Denote(22)
Δ
f
j
=
f
(
x
j
)
-
f
(
x
j
+
d
j
)
,
(23)
Δ
q
j
=
q
j
(
0
)
-
q
j
(
d
j
)
,
(24)
r
j
=
Δ
f
j
Δ
q
j
;
the general trust region algorithm [17] for solving the unconstrained optimization problem (20) is given as follows.Algorithm 1.
Step
0. Given a starting point x
0
∈
R
n, Δ
0
>
0 is the initial trust region radium, ɛ
>
0, B
0
=
Hesse
(
x
0
). Set j
:
=
0.
Step
1. Calculate g
j, if g
j
≤
ɛ, and stop iteration; otherwise, go to Step 2.
Step
2. Utilize the smoothing Newton method to obtain the solution d
j of the subproblem (21).
Step
3. Calculate the r
j of formula (24).
Step
4. Regulate trust region radius. If r
j
<
0.25, let Δ
j
+
1
:
=
0.5
Δ
j; if r
j
>
0.75 and d
j
=
Δ
j, let Δ
j
+
1
:
=
2
Δ
j; otherwise, let Δ
j
+
1
:
=
Δ
j.
Step
5. If r
j
>
0.25, let x
j
+
1
:
=
x
j
+
d
j, update B
j to be B
j
+
1
=
Hesse
(
x
j
+
1
), let j
:
=
j
+
1, and go to Step 1; otherwise, let x
j
+
1
=
x
j, j
:
=
j
+
1, and go to Step 2, where Δ
0
=
1
/
10
g
(
x
0
), ɛ
=
1
0
-
6.Obviously, the above trust region subproblem (21) is an unconstrained optimization problem whose objective function is a quadratic. For the sake of convenience, we first denote
(25)
z
=
μ
,
λ
,
d
,
(26)
Y
(
z
)
:
=
μ
λ
-
d
2
2
+
Δ
j
2
-
λ
+
d
2
2
-
Δ
j
2
2
+
4
μ
2
B
j
+
λ
E
d
-
g
j
,
(27)
β
(
z
)
=
γ
Y
z
min
⁡
1
,
Y
z
,
respectively. Then, the following smoothing Newton algorithm [17] is given to solve the subproblem (21).Algorithm 2.
Step
0. Given d
0
∈
R
n, z
0
=
(
μ
0
,
λ
0
,
d
0
), z
~
0
=
(
μ
0
,
0,0
). Select parameters δ
,
σ
,
γ
∈
(
0,1
), μ
0
γ
<
1, and γ
Y
(
z
0
)
<
1; set h
:
=
0.
Step
1. Calculate Y
(
z
h
), if Y
(
z
h
)
=
0, and terminate algorithm; otherwise, calculate β
h
=
β
(
z
h
).
Step
2. Obtain the solution for equations Y
(
z
h
)
+
Y
′
(
z
h
)
Δ
z
h
=
β
h
z
~, and the solution is Δ
z
h
=
(
Δ
μ
h
,
Δ
λ
h
,
Δ
d
h
).
Step
3. Suppose that m
h is the minimum nonnegative integer meeting Y
(
z
h
+
δ
m
h
Δ
z
h
)
≤
[
1
-
σ
(
1
-
β
μ
0
)
δ
m
h
]
Y
(
z
h
). Let α
h
:
=
δ
m
h and z
h
+
1
=
z
h
+
α
h
Δ
z
h.
Step
4. Let h
:
=
h
+
1, and go to Step 1.


## 3.2. A New Newton-Type Trust Region Algorithm

Optimizing ability of Newton-type trust region algorithm (N-TTRA) highly depends on initial parameters. And the algorithm is short of global optimizing capability although its local optimizing is fast. In order to enhance the global optimization capability and get rid of dependence on initial parameters, we construct a new Newton-type trust region algorithm by combining a new BFGS updating formula with N-TTRA. The basic idea is to updateB
j
+
1 with new BFGS formulas at Step 5 in Algorithm 1. This idea guarantees positive definiteness of B
j
+
1 to improve the global optimization capability; meanwhile, it gets rid of dependence on the initial parameter selection. New adjustment formula of BFGS is as follows:
(28)
B
j
+
1
=
B
j
+
y
j
*
y
j
*
T
s
j
T
y
j
*
-
B
j
s
j
s
j
T
B
j
s
j
T
B
j
s
j
,
where y
j
*
=
y
j
T
s
j
/
y
j
T
s
j
y
j, s
j
=
x
j
+
1
-
x
j, and y
j
=
g
j
+
1
-
g
j. Here, y
j
* guarantees positive definiteness of B
j
+
1, so the new Newton-type trust region algorithm (NN-TTRA) uses (28) as the new adjustment formula of B
j
+
1 at Step 5 in Algorithm 1. It is proved that the method is of super-linear convergence [22]. The flow chart of NN-TTRA is given in Figure 1.Figure 1
Flow chart of NN-TTRA.


## 4. Numerical Experiments and Results Analysis

In this paper, we provide actual observed data of flood runoff process between Chenggouwan and Linqing segment in Nanyunhe River of Haihe River Basin. (Length of the reach is 83.8 km, where there is no tributary, but a levee control on both sides. There may occur lifting irrigation during the water delivery, and flood water may discharge into the reach when rainfall is high. But these situations have little effect on flood, where the routing time intervalΔ
t
=
12 h.) The detailed data can be seen in [23]. Here, we will give the numerical experiments from the following three aspects.Firstly, in order to verify the advantage of our new parameter estimation models (16) by using the generalized trapezoid formula to approximate (1), we use the NN-TTRA to solve the above unconstrained optimization problem (16) by choosing different θ. Furthermore, we get the corresponding parameters K and x for each θ. From these parameters, we use (14) to obtain the calculated outflow Q
~
(
t
i
), where Q
~
(
t
1
)
=
Q
(
t
1
). At last, the average absolute errors (AAE) and the average relative errors (ARE) can be given as follows:
(29)
AAE
=
1
n
∑
i
=
1
n
‍
Q
~
(
t
i
)
-
Q
(
t
i
)
,
ARE
=
1
n
∑
i
=
1
n
‍
Q
~
(
t
i
)
-
Q
(
t
i
)
Q
(
t
i
)
.
From (29), Tables 1, 2, and 3 list the AAE and ARE for different θ for flood routing in 1960, 1961, and 1964, respectively. It is shown from Tables 1–3 that the AAE and ARE are the smallest when θ
=
1
/
3. In other words, the numerical results given in Tables 1–3 confirm the theoretical analysis presented in Section 2.2.Table 1
Numerical results with differentθ for flood routing in 1960.
θ
AAE
ARE (%)
θ
AAE
ARE (%)
θ
=
0
8.208160
2.462451
θ
=
1
/
3
8.208122
2.462438
θ
=
1
/
10
8.208159
2.462451
θ
=
1
/
2
8.208163
2.462452
θ
=
1
/
7
8.208154
2.462449
θ
=
2
/
3
8.208162
2.462452
θ
=
1
/
5
8.208169
2.462454
θ
=
1
8.208162
2.462451Table 2
Numerical results with differentθ for flood routing in 1961.
θ
AAE
ARE (%)
θ
AAE
ARE (%)
θ
=
0
4.011335
0.998421
θ
=
1
/
3
4.011335
0.998421
θ
=
1
/
10
4.011335
0.998421
θ
=
1
/
2
4.011336
0.998421
θ
=
1
/
7
4.011336
0.998421
θ
=
2
/
3
4.011336
0.998421
θ
=
1
/
5
4.011336
0.998421
θ
=
1
4.011335
0.998421Table 3
Numerical results with differentθ for flood routing in 1964.
θ
AAE
ARE (%)
θ
AAE
ARE (%)
θ
=
0
10.949204
2.662713
θ
=
1
/
3
10.949203
2.662712
θ
=
1
/
10
10.949205
2.662713
θ
=
1
/
2
10.949204
2.662713
θ
=
1
/
7
10.949208
2.662713
θ
=
2
/
3
10.949203
2.662712
θ
=
1
/
5
10.949205
2.662713
θ
=
1
10.949204
2.662712Secondly, to illustrate the efficiency of the NN-TTRA presented in this paper, for different parameter initial values, we employ the N-TTRA and the NN-TTRA to solve the above unconstrained optimization problem (16), respectively, where the parameter θ is set to 1
/
3 and a maximum number of iterations of the N-TTRA and the NN-TTRA are set to 150. The numerical results are shown in Table 4. From Table 4, we can indicate that the NN-TTRA not only gets rid of dependence on the initial values of parameters, but also guarantees that the number of iterations is finite. Meanwhile, through this algorithm the global optimum value is obtained. However, the N-TTRA highly depends on initial values of parameters. When the selected initial value is far from the optimal ones, the algorithm hardly finds the global optimum and even cannot complete iterations in finite times.Table 4
Comparison of performance of the two algorithms.
Algorithm
Initial value
Iterations
K
x
N-TTRA
(1, 1)
5
0
−1.938700
(2, 2)
—
—
—
(11, 1)
3
11.190740
0.998456
NN-TTRA
(1, 1)
27
11.190740
0.998456
(2, 2)
25
11.190740
0.998456
(11, 1)
4
11.190740
0.998456Finally, Table5 lists the numerical results calculated by using NN-TTRA to solve problem (16), where θ
=
1
/
3. Meanwhile Table 5 also gives the numerical results obtained by using method of trial-and-error (TAE) [2], the least-square method (L-SM) [6], and direct optimal method (DOM) [18], respectively. Here, the observed data is the outflow of flood runoff between Chenggouwan and Linqing segment in Nanyunhe River of Haihe Basin in 1960. It can be seen from Table 5 that the AAE and ARE of NN-TTRA are all less than those of other estimation methods. Thus, we conclude that it is very effective by using NN-TTRA to estimate the parameters of Muskingum model. In addition, the calculated values and observed values of the flows in 1960, 1961, and 1964 are plotted in Figures 2, 3, and 4, respectively. From Figures 2–4, it can be seen that the calculated flow data via the NN-TTRA highly coincides with the observed flow data. In brief, the accuracy of the method is satisfactory.Table 5
Comparison of results via several parameter estimation methods.
Algorithm
K
x
AAE
ARE (%)
TAE [2]
12.400000
0.100000
10.230000
3.150000
L-SM [6]
11.790000
−0.325000
9.800000
2.960000
DOM [18]
12.440000
−0.260000
9.700000
2.920000
NN-TTRA
11.190740
0.998456
8.208122
2.462438Figure 2
Comparison between calculated values and observed values in 1960.Figure 3
Comparison between calculated values and observed values in 1961.Figure 4
Comparison between calculated values and observed values in 1964.


## 5. Conclusions

By combining a new BFGS adjustment formula with N-TTRA, this paper implements parameter estimation of Muskingum model by selecting different initial values for parameters and comparing it with N-TTRA. The results show that the NN-TTRA can get rid of the dependence on initial value selection for parameters when searching solutions for Muskingum model parameters, which avoids the influence of initial value selection for parameters on the optimization results and averts local optimum. What is more, this algorithm is of application value in flood disaster management and should be generalized. In addition, this algorithm can be extended to other similar parameter estimation problems to help obtain excellent results.



---

*Source: 634852-2014-12-21.xml*