# An Improved Global Harmony Search Algorithm for the Identification of Nonlinear Discrete-Time Systems Based on Volterra Filter Modeling

Authors: Zongyan Li; Deliang Li
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

This paper describes an improved global harmony search (IGHS) algorithm for identifying the nonlinear discrete-time systems based on second-order Volterra model. The IGHS is an improved version of the novel global harmony search (NGHS) algorithm, and it makes two significant improvements on the NGHS. First, the genetic mutation operation is modified by combining normal distribution and Cauchy distribution, which enables the IGHS to fully explore and exploit the solution space. Second, an opposition-based learning (OBL) is introduced and modified to improve the quality of harmony vectors. The IGHS algorithm is implemented on two numerical examples, and they are nonlinear discrete-time rational system and the real heat exchanger, respectively. The results of the IGHS are compared with those of the other three methods, and it has been verified to be more effective than the other three methods on solving the above two problems with different input signals and system memory sizes.


---

## Body

## 1. Introduction

The Volterra model is a kind of nonlinear filter model, which is usually employed to track and identify plenty of complex nonlinear systems. In order to enhance the quality of the system identification, it is very crucial issue to select optimum model coefficient called the kernel. Therefore, the Volterra model is essentially an extension of linear filter model to nonlinear case. During the past decade, there have been many research works on the Volterra model. Campello et al. [1] tackled the problem of expanding Volterra models using Laguerre functions. The global optimal solution is obtained when each multidimensional kernel of the model is decomposed into a number of independent orthonormal bases. Furthermore, the solution obtained is able to minimize the upper bound of the squared norm of the error resulting from the practical truncation of the Laguerre series expansion into a finite number of functions. Masugi and Takuma [2] described a Volterra system-based nonlinear study of video-packet transmission over IP networks. Based on the Volterra system, the authors performed a time-series analysis of measured data for network response evaluation. The novel method can reproduce the time-series responses observed in video-packet transmission over the Internet, characterizing nonlinear dynamic behaviors such that the obtained results gave an appropriate depiction of network conditions at different times. Gruber et al. [3] presented a nonlinear model predictive control (NMPC) method based on a second-order Volterra series model for greenhouse temperature control using natural ventilation. These models, denoting the simple and logical extension of convolution models, are capable of describing the nonlinear dynamic feature of the ventilation and other environmental conditions on the greenhouse temperature. Many applications of Volterra series modeling were executed in the frequency-domain based on the Generalized Frequency Response Functions (GFRF) [4, 5]. In the light of the wide applications of Volterra series, Li and Billings [6] presented an approach to estimate the GFRF in a piecewise manner for duffing type oscillators for the underrepresented weakly nonlinear region. Additionally, they devised a new method to obtain the energy contributions of each order kernel. The new nonparametric method can not only construct the amplitude-invariant GFRF over a certain excitation range, but also avoid building large sets of time-domain models. In addition, the Volterra filter model can be found in some other application areas [7–14].Chang [7] devised an improved particle swarm optimization (IPSO) to implement system identification based on Volterra filter model. In this paper, we develop an improved global harmony search (IGHS) algorithm and try the IGHS as an efficient candidate for system identification based on Volterra filter model. The harmony search (HS) algorithm was firstly proposed by [15]. The HS is a simple but efficient algorithm, and its many improved versions have been applied into many problems including reliability problems [16], reactor core fuel management optimization [17], and sizing optimization of truss structures [18].The paper is organized as follows. In Section2, the pruned second-order Volterra filter model is simply presented. In Section 3, a novel global harmony search algorithm is introduced. In Section 4, an improved global harmony search algorithm is proposed, and its procedure is fully explained. In Section 5, four harmony search algorithms are used for two examples with different signal inputs and memory sizes. We end this paper with some conclusions in Section 6.


## 2. Volterra Filter Model and Its Pruned Second-Order Form

The Volterra filter model is an efficient method for the identification of nonlinear discrete systems, and it has come into researchers’ notice in recent decades. The discrete form of Volterra filter model of theqth order [7] is given by(1)
y
n
=
h
0
+
∑
k
1
=
0
N
-
1
h
1
k
1
x
n
-
k
1
+
∑
k
1
=
0
N
-
1
∑
k
2
=
0
N
-
1
h
2
k
1
,
k
2
x
n
-
k
1
x
n
-
k
2
+
⋯
+
∑
k
1
=
0
N
-
1
∑
k
2
=
0
N
-
1
⋯
∑
k
q
=
0
N
-
1
h
q
k
1
,
k
2
,
…
,
k
q
∏
i
=
1
q
x
n
-
k
i
=
h
0
+
∑
j
=
1
q
∑
k
1
=
0
N
-
1
⋯
∑
k
j
=
0
N
-
1
h
j
k
1
,
k
2
,
…
,
k
j
∏
i
=
1
j
x
n
-
k
i
,where N denotes the system memory size. Equation (1) denotes the Volterra filter model with the infinite series. However, this model is hard to compute and master due to its complex and expatiatory formula. In this paper, we only study its simplified and approximate form called the truncated second-order Volterra model [7, 19], which is stated as follows:(2)
y
n
=
h
0
+
∑
k
=
1
N
h
k
x
n
-
k
+
1
+
∑
k
1
=
0
N
-
1
∑
k
2
=
k
1
N
-
1
h
k
1
,
k
2
x
n
-
k
1
x
n
-
k
2
.In order to facilitate expression, (2) can also be expressed as the following vector form:(3)
y
n
=
H
X
T
.Here, the superscriptT represents the transpose of a vector and H stands for the Volterra kernel vector given by(4)
H
=
h
0
,
h
1
,
…
,
h
N
,
h
0,0
,
h
0,1
,
…
,
h
0
,
N
-
1
,
h
1,1
,
…
,
h
N
-
1
,
N
-
1
.In addition,X denotes the Volterra input vector given by(5)
X
=
1
,
x
n
,
…
,
x
n
-
N
+
1
,
x
2
n
,
x
n
x
n
-
1
,
…
,
x
n
x
n
-
N
-
1
,
x
2
n
-
1
,
…
,
x
2
n
-
N
-
1
.In the light of (3), the vector lengths of both H and X are the same and are calculated as follows [19]:(6)
L
=
C
N
+
2
N
=
1
+
N
+
N
N
+
1
2
=
N
+
1
N
+
2
2
.To achieve the nearest approximation of the actual system output, appropriate kernel vectorH should be determined under the input vector X. In this paper, an improved global harmony search (IGHS) algorithm is proposed to determine kernel vector H. The IGHS is an improved version of novel global harmony search (NGHS) algorithm [20]; thus, both the NGHS and the IGHS will be presented in the following sections.


## 3. Novel Global Harmony Search (NGHS) Algorithm

Novel global harmony search (NGHS) algorithm [20] is a variant of harmony search (HS) algorithm [15], and it is superior to the HS for solving unconstrained optimization problems. The NGHS improvises new harmony vectors by combining position updating and mutation. Concretely, the steps of NGHS are explained as follows.Step 1 (initialize the NGHS parameters and the problem parameters).
The NGHS parameters consist of harmony memory sizeH
M
S, the number of improvisations N
I, and mutation rate p
m. In addition, the problem parameters include the number of problem variables L, the lower bound x
_
j, and upper bound x
¯
j of the jth (j
=
1,2
,
…
,
L) problem variable. Furthermore, the number of improvisations (N
I) is actually the total number of generations for adjusting the parameters related to the identification problem. In each generation, only one new candidate solution including L parameters is generated, and this solution is accepted if and only if it is better than the worst one of the previous solutions.Step 2 (initialize harmony memory (HM)).
The initial harmony memory (H
M) can be expressed in the following matrix form:(7)
H
M
=
x
1
1
x
2
1
⋯
x
L
1
x
1
2
x
2
2
⋯
x
L
2
⋮
⋮
⋱
⋮
x
1
H
M
S
x
2
H
M
S
⋯
x
L
H
M
S
,where x
j
i stands for the jth (j
=
1
,
…
,
L) variable of the ith (i
=
1
,
…
,
H
M
S) harmony vector. Moreover, it is randomly produced from a uniform distribution in the ranges [
x
_
j
,
x
¯
j
]
(
j
=
1,2
,
…
,
L
). In HM, any vector x
i
=
(
x
1
i
,
x
2
i
,
…
,
x
L
i
) (i
=
1
,
…
,
H
M
S) represents a candidate solution of the parameters (as in (4)) needed for solving the identification problem. Furthermore, the length of x
i is equal to L, which is exactly the number of the parameters in the Volterra kernel vector H (as in (6)).Step 3 (improvise a new harmony).
Improvisation is actually the operation of producing a new harmony vector. For the NGHS, its improvisation mainly includes two steps, and they are position updating and genetic mutation, respectively. More specifically, the improvisation can be presented in Table1.
Here, “best" and “worst” stand for the indexes of the best harmony and the worst harmony inH
M, respectively. r
1, r
2, and r
3 denote three uniformly generated random numbers in [
0,1
]. With respect to position updating, a new harmony vector is improvised near the best harmony vector, which can facilitate the convergence rate of the NGHS. On the other hand, it is worth noticing that genetic mutation is an event of small probability, and it is utilized to avoid the premature convergence of the NGHS.Table 1
The computational procedure of the NGHS improvisation.
Line
The improvisation of the NGHS
(1)
For
j
=
1 to L
(2)
x
R
=
2
×
x
j
b
e
s
t
-
x
j
w
o
r
s
t
(3)
If
x
R
>
x
¯
j
(4)
x
R
=
x
¯
j
(5)
ElseIf
x
R
<
x
_
j
(6)
x
R
=
x
_
j
(7)
EndIf
(8)
x
j
′
=
x
j
w
o
r
s
t
+
r
1
×
(
x
R
-
x
j
w
o
r
s
t
) % Position updating
(9)
If
r
2
≤
p
m
(10)
x
j
′
=
x
_
j
+
r
3
×
(
x
¯
j
-
x
_
j
) % Genetic mutation
(11)
EndIf
(12)
EndForStep 4 (update harmony memory).
Replace the worst harmony vectorx
w
o
r
s
t of H
M with the new improvised harmony vector x
′ no matter whether x
′ is better or worse than x
w
o
r
s
t.Step 5 (check the stopping criterion).
Steps3 and 4 are repeated until the number of improvisations (N
I) is reached.


## 4. An Improved Global Harmony Search (IGHS) Algorithm

In this paper, an improved global harmony search (IGHS) algorithm is proposed to implement the system identification based on Volterra filter model. The IGHS is an improved version of the NGHS, and it is different from the NGHS in the following two aspects.(
1) The Modification of Genetic Mutation. The genetic mutation of the NGHS is conducted according to uniform distribution. To further explore and exploit the solution space, the uniform distribution is replaced with the other two probability distributions, and they are normal distribution [21–23] and Cauchy distribution [22, 23], respectively. Therefore, the new genetic mutation are implemented by using the following two equations:(8)
x
j
′
=
r
a
n
d
n
j
x
j
b
e
s
t
,
0.1
,
x
j
′
=
r
a
n
d
c
j
x
j
b
e
s
t
,
0.1
.Here,j denotes the index of the jth component of the improvised harmony vector, and x
j
b
e
s
t denotes the jth component of the best harmony vector in H
M. r
a
n
d
n
j
(
x
j
b
e
s
t
,
0.1
) stands for the normal distribution with mean x
j
best and standard deviation 0.1, and r
a
n
d
c
j
(
x
j
b
e
s
t
,
0.1
) stands for the Cauchy distribution with location parameter x
j
b
e
s
t and scale parameter 0.1. If overflow happens to x
j
′ generated using either normal distribution or Cauchy distribution, it will be truncated to [
x
¯
j
,
x
_
j
]. On the other hand, if the condition r
2
<
p
m (as in Table 1) is satisfied, either normal distribution or Cauchy distribution will be utilized to carry out genetic mutation, and the probability of using any one of the two distributions is equal to 0.5. By adopting normal distribution and Cauchy distribution, the capacity of escaping from the local optimums is enhanced for the IGHS. In the meantime, the solution space can be fully explored and exploited, which is beneficial for improving the quality of the improvised harmony vector x
′.In addition to the equation of genetic mutation, we also modify the mutation rate (p
m) of the NGHS. The parameter p
m is used to determine whether or not a variable of x
′ adopts genetic mutation. A large p
m value is beneficial for searching in a large scope, while a small p
m value is helpful to search in a small scope. To well balance the global search and local search of the IGHS, the p
m value decreases dynamically with increasing generations (N
I) as follows:(9)
p
m
t
=
p
¯
m
-
p
¯
m
-
p
_
m
N
I
0
×
t
,
if
t
≤
N
I
0
;
p
_
m
,
if
t
>
N
I
0
, where p
_
m and p
¯
m represent the minimal and maximal mutation rates. t (1
≤
t
≤
N
I) denotes the current generation number, and the IGHS is stopped when the current generation number t reaches the maximal generation number N
I. In addition, N
I
0 denotes a fixed generation number, and it is set to 3
N
I
/
4 here.(
2) Introduce and Modify an Opposition-Based Learning (OBL) [24]. In order to improve the convergence of the IGHS, a method called opposition-based learning (OBL) [24] is firstly introduced. In short, this technique can be stated as follows:(10)
o
x
j
i
=
x
j
L
+
x
j
U
-
x
j
i
.Here,o
x
i is the opposite solution of x
i, and o
x
j
i is the jth variable of o
x
i. Furthermore, x
j
L and x
j
U represent the minimal and the maximal values in the set {
x
j
1
,
x
j
2
,
…
,
x
j
H
M
S
}, and they are expressed as follows:(11)
x
j
L
=
m
i
n
1
≤
i
≤
H
M
S
x
j
i
,
x
j
U
=
m
a
x
1
≤
i
≤
H
M
S
x
j
i
.According to (10), only the opposite value of x
j
i is employed to produce the corresponding variable in the range of [
x
j
L
,
x
j
U
]. In other words, once the x
j
i value is determined, it has an exclusive opposite value o
x
j
i which is equal to x
j
L
+
x
j
U
-
x
j
i.In this paper, we enhance the IGHS improvisation by modifying the OBL technique. More specifically, theo
x
j
i value is randomly produced in the range of [
x
j
L
,
x
j
U
], and it is given by(12)
o
x
j
i
=
x
j
L
+
r
a
n
d
×
x
j
U
-
x
j
L
,where rand
(
) denotes a random number uniformly generated in [
0,1
]. By using (12), the IGHS can yield numerous possible values in the range of [
x
j
L
,
x
j
U
]. Compared to the OBL technique, its modified version has wider searching space, which is in favor of fine-tuning of the best harmony vectors.Based on (10)-(11), a new harmony memory H
M
′ can be generated, and it is given by(13)
H
M
′
=
o
x
1
1
o
x
2
1
⋯
o
x
L
1
o
x
1
2
o
x
2
2
⋯
o
x
L
2
⋮
⋮
⋱
⋮
o
x
1
H
M
S
o
x
2
H
M
S
⋯
o
x
L
H
M
S
.Here,o
x
i
=
(
o
x
1
i
,
o
x
2
i
,
…
,
o
x
L
i
) stands for the ith (i
=
1
,
…
,
H
M
S) harmony vector in H
M
′, and o
x
j
i denotes its jth (j
=
1
,
…
,
L) variable. After performing the modified OBL operation, each o
x
i is compared with the corresponding harmony vector x
i. More precisely, if o
x
i is better than x
i, x
i should be replaced with o
x
i. It is worth noticing that the modified OBL is an event of small probability. In other words, this technique is an auxiliary step of the IGHS improvisation, and it is mainly used to improve the quality of the improvised harmony vector.Based on the above detailed illustrations about the two modifications of the NGHS improvisation, the IGHS improvisation can be summarized in Table2.Table 2
The computational procedure of the IGHS improvisation.
Line
The improvisation of the IGHS
(1)
For
j
=
1 to L
(2)
x
R
=
2
×
x
j
b
e
s
t
-
x
j
w
o
r
s
t
(3)
If
x
R
>
x
¯
j
(4)
x
R
=
x
¯
j
(5)
ElseIf
x
R
<
x
_
j
(6)
x
R
=
x
_
j
(7)
EndIf
(8)
x
j
′
=
x
j
w
o
r
s
t
+
r
1
×
(
x
R
-
x
j
w
o
r
s
t
)% Position updating
(9)
If
r
2
≤
p
m
t%  p
m
t is calculated according to (9)
(10)
If
r
3
≤
0.5
(11)
x
j
′
=
rand
⁡
n
j
(
x
j
b
e
s
t
,
0.1
)% Normal distribution
(12)
Else
(13)
x
j
′
=
rand
⁡
c
j
(
x
j
b
e
s
t
,
0.1
)% Cauchy distribution
(14)
EndIf
(15)
If
x
j
′
>
x
¯
j
(16)
x
j
′
=
x
¯
j
(17)
ElseIf
x
j
′
<
x
_
j
(18)
x
j
′
=
x
_
j
(19)
EndIf
(20)
EndIf
(21)
EndFor
(22)
If
r
4
<
0.1%  Modified OBL operation with a probability of 0.1
(23)
For
i
=
1 to H
M
S
(24)
For
j
=
1 to L
(25)
o
x
j
i
=
x
j
L
+
rand
⁡
(
)
×
(
x
j
U
-
x
j
L
)
(26)
EndFor
(27)
If
f
(
o
x
i
)
<
f
(
x
i
)
(28)
x
i
=
o
x
i
(29)
Endif
(30)
EndFor
(31)
EndIfBy utilizing the IGHS, the identification diagram of nonlinear discrete-time systems based on the pruned second-order Volterra model is shown in Figure1. Besides, the nomenclatures appearing in this diagram are explained in Table 3.Table 3
Nomenclatures used in Volterra filter modeling of nonlinear discrete-time system.
x
[
n
]
The digital input signal
y
[
n
]
The output signal of the second-order Volterra filter model (as (2))
d
[
n
]
The output of the unknown nonlinear discrete-time system
r
[
n
]
The measurement noise
y
^
[
n
]
The sum ofd
[
n
] and r
[
n
]
e
[
n
]
The error signal betweeny
^
[
n
] and y
[
n
]Figure 1
Volterra filter modeling of nonlinear discrete-time system using the IGHS algorithm.The goal of the IGHS is to find the optimal kernel vectorH for Volterra filter model so that the difference between the estimated output y
[
n
] and the actual system output y
^
[
n
] is minimized. Thus, it is advisable to find an objective function to satisfy the design requirement. Here we adopt the following objective function [7]:(14)
E
e
2
n
=
1
T
∑
n
=
0
T
-
1
e
2
n
=
1
T
∑
n
=
0
T
-
1
y
^
n
-
y
n
2
.Here, E
[
e
2
[
n
]
] stands for the mean square error (MSE), and T denotes the sampling number. By using the IGHS to optimize MSE, we can obtain the most appropriate kernel vector for the second-order Volterra model.


## 5. Experimental Results and Analysis

To verify the validity of the IGHS on identifying nonlinear system based on second-order Volterra filter model, two examples including the highly nonlinear discrete-time rational system and the real heat exchanger are considered. More specifically, these two examples are explained as follows.


### 5.1. Example 1

The first example is the highly nonlinear discrete-time rational system [7], and its mathematical model is stated as follows:(15)
d
n
=
0.3
d
2
n
-
1
+
0.8
x
n
-
1
+
0.6
d
n
-
2
1
+
x
2
n
-
1
+
d
2
n
-
1
.Here, two types of input signals [7] are considered. The first is defined as Example 1a whose input is a random signal uniformly produced in [
-
1,1
]. The second is defined as Example 1b whose input is a cosine signal x
[
n
]
=
0.8
c
o
s
(
(
π
/
9
)
n
). For both Examples 1a and 1b, the measurement noise r
[
n
] is always supposed to be a Gaussian noise of N
(
0,0.001
).In this experiment, the IGHS is used for Example1a with N
=
5. Moreover, the IGHS parameters are set as follows: the maximal mutation rate p
¯
m
=
0.1, the minimal mutation rate p
_
m
=
0.1, harmony memory size H
M
S
=
5, the number of improvisations N
I
=
2000, and sampling number T
=
100. Moreover, T is actually the total number of output samples (actual system output or estimated output). According to (14), there are T sampling values for actual system output y
^
[
n
] (1
≤
n
≤
T). After implementing the IGHS, we can obtain the comparison between the estimated output y
[
n
] and the actual output y
^
[
n
] in Figure 2.Figure 2
Comparison of actual system output and Volterra model output for Example1a (N
=
5).From Figure2, it can be seen that the estimated output y
[
n
] approximates the actual output y
^
[
n
] well. In addition, the minimal mean square error (MSE) obtained by the IGHS is equal to 5.4594
e
-
003, which is a satisfactory result. Besides, we use Example 1a with N
=
8 to investigate the performance of the IGHS for second-order Volterra filter model with large memory size. The parameters of the IGHS are the same as those used for Example 1a (N
=
5) except N
I, and it is set to 5000 in this experiment. By performing the IGHS, we can obtain the comparison between the estimated output y
[
n
] and the actual output y
^
[
n
] in Figure 3.Figure 3
Comparison of actual system output and Volterra model output for Example1a (N
=
8).It is clear from Figure3 that a satisfactory approximation can be attained by utilizing the IGHS. Additionally, the minimal MSE yielded by the IGHS is equal to 3.7074
e
-
003 for Example 1a with N
=
8, which provides better modeling capacity.In addition to random signal, another testing input signalx
[
n
]
=
0.8
c
o
s
(
(
π
/
9
)
n
) is used to investigate the performance of second-order Volterra filter model using the IGHS. For Example 1b with N
=
5 and N
=
8, the IGHS parameters are the same as those for Example 1a, and Figures 4 and 5 display the comparisons of results for Example 1b with N
=
5 and N
=
8. As expected, the difference between the estimated output y
[
n
] and the actual output y
^
[
n
] is small in each case. Additionally, two satisfying MSEs can be obtained, respectively, for N
=
5 and N
=
8 by carrying out the IGHS, and they are equal to 1.71809876
e
-
003 and 1.7119
e
-
003, respectively.Figure 4
Comparison of actual system output and Volterra model output for Example1b (N
=
5).Figure 5
Comparison of actual system output and Volterra model output for Example1b (N
=
8).


### 5.2. Example 2

The second example is the real heat exchanger and its mathematical model [7, 25] is given by(16)
w
n
=
x
n
-
1.3228
x
2
n
+
0.7671
x
3
n
-
2.1755
x
4
n
,
(17)
d
n
=
-
6.5306
z
-
1
+
5.5652
z
-
2
1
-
1.608
z
-
1
+
0.6385
z
-
2
w
n
, where x
[
n
] is the process input denoting the flow rate and is constrained by the range of [
0,1
], w
[
n
] stands for the static nonlinearity, and d
[
n
] stands for the process output temperature. By combining (16) and (17), the simplified mathematical model can be expressed as the following difference equation:(18)
d
n
=
1.608
d
n
-
1
-
0.6385
d
n
-
2
-
6.5306
w
n
-
1
+
5.5652
w
n
-
2
.Here, two types of input signals [7, 25] are considered. The first is defined as Example 2a whose input is randomly generated in the range [
0.1,0.9
]. The second is defined as Example 2b whose input is a sine signal x
[
n
]
=
0.4
s
i
n
(
(
π
/
6
)
n
)
+
0.5. According to [7, 25], the measurement noise is excluded from the problem model; thus, we have r
[
n
]
=
0.In this experiment, the memory sizeN is set to 8 for second-order Volterra filter model. To fulfill the physical input requirement, the modeling input x
[
n
] is randomly generated in the range [
0.1,0.9
] and the testing input is then set to x
[
n
]
=
0.4
s
i
n
(
(
π
/
6
)
n
)
+
0.5. Moreover, the IGHS parameters used for Example 2 are the same as those of Example 1. By implementing the IGHS, the outputs of second-order Volterra model are shown in Figures 6 and 7, respectively, for the modeling input and testing input. Based on the careful observations on experimental results, we can confirm that the satisfactory approximation is attained in each case.Figure 6
Comparison of actual system output and Volterra model output for Example2a (N
=
8).Figure 7
Comparison of actual system output and Volterra model output for Example2b (N
=
8).


### 5.3. Comparison of the IGHS with the Other Three HSs

In this paper, four methods consisting of the HS [15], IHS [26], NGHS [20], and IGHS are used to solve the above problems with different input signals and system memory sizes. The parameters of these four methods are set as follows. For the HS, harmony memory considering rate H
M
C
R
=
0.95, pitch adjusting rate P
A
R
=
0.3, and bandwidth b
w
=
0.01. For the IHS, H
M
C
R
=
0.95, the minimal pitch adjusting rate P
A
R
m
i
n
=
0.35, the maximal pitch adjusting rate P
A
R
m
a
x
=
0.99, the minimal bandwidth b
w
m
i
n
=
1
0
-
5, and the maximal bandwidth b
w
m
a
x
=
0.1. For the NGHS, mutation rate p
m
=
5
×
1
0
-
2. For the IGHS, the minimal mutation rate p
_
m
=
0.01, the maximal mutation rate p
¯
m
=
0.1, and N
I
0
=
3
N
I
/
4. Note that all the methods are executed under identical conditions (harmony memory size and the number of improvisations). More specifically, the harmony memory size (H
M
S) of each method is set to 5, and the N
I (Number of improvisations) value of each method is set to 2000 and 5000, respectively, for the system memory sizes N
=
5 and N
=
8. With respect to second-order Volterra filter model, sampling number is set to T
=
100, and the range of each variable of the kernel vector is set as [
-
1,1
]. Matlab 7.0 is used to execute the above procedure under the environment of Intel(R) Core(TM) i5-2410M CPU @ 2.30 GHz. 20 independent runs are performed in each case, and the optimization results are presented in Table 4.Table 4
Comparison of HS, IHS, NGHS, and IGHS on six problems.
Problem
NI
Algorithm
ACT
Best
Worst
Mean
Std
Example1a (N
=
5)
2000
HS
1.659
1.1722
e
-
002
3.5945
e
-
002
1.9427
e
-
002
8.0506
e
-
003
IHS
1.7011
8.4989
e
-
003
2.6525
e
-
002
1.6498
e
-
002
5.9693
e
-
003
NGHS
1.781
5.9558
e
-
003
3.0649
e
-
002
1.2990
e
-
002
7.8612
e
-
003
IGHS
2.6365
5.4594
e
-
003
7.9210
e
-
003
5.9763
e
-
003
6.6092
e
-
004
Example1a (N
=
8)
5000
HS
4.3975
2.9721
e
-
002
9.3845
e
-
002
5.4468
e
-
002
1.5500
e
-
002
IHS
4.4121
3.1554
e
-
002
8.1672
e
-
002
5.1982
e
-
002
1.4302
e
-
002
NGHS
4.4334
5.9619
e
-
003
4.9489
e
-
002
1.5350
e
-
002
9.0064
e
-
003
IGHS
6.6589
3.7074
e
-
003
1.0237
e
-
002
5.6536
e
-
003
1.7608
e
-
003
Example1b (N
=
5)
2000
HS
0.57743
1.72096501
e
-
003
4.06941185
e
-
003
2.15386021
e
-
003
5.99250994
e
-
004
IHS
0.56219
1.71809880
e
-
003
1.82397928
e
-
003
1.72847844
e
-
003
2.71096253
e
-
005
NGHS
0.56886
1.71809919
e
-
003
1.83200101
e
-
003
1.72490455
e
-
003
2.54479841
e
-
005
IGHS
0.84918
1.71809876
e
-
003
1.71833108
e
-
003
1.71811926
e
-
003
5.18365165
e
-
008
Example1b (N
=
8)
5000
HS
2.6544
1.7202
e
-
003
2.3210
e
-
003
1.8309
e
-
003
1.5559
e
-
004
IHS
2.651
1.7119
e
-
003
1.7721
e
-
003
1.7194
e
-
003
1.8286
e
-
005
NGHS
2.6787
1.7119
e
-
003
1.8486
e
-
003
1.7188
e
-
003
3.0568
e
-
005
IGHS
4.0572
1.7119
e
-
003
1.7119
e
-
003
1.7119
e
-
003
1.4972
e
-
009
Example2a (N
=
8)
5000
HS
5.3252
3.8006
e
+
000
4.0102
e
+
000
3.8879
e
+
000
5.2382
e
-
002
IHS
5.2646
3.7901
e
+
000
3.9328
e
+
000
3.8548
e
+
000
4.3294
e
-
002
NGHS
5.3327
3.6604
e
+
000
3.7327
e
+
000
3.6760
e
+
000
1.6444
e
-
002
IGHS
7.9408
3.6560
e
+
000
3.6628
e
+
000
3.6581
e
+
000
1.9936
e
-
003
Example2b (N
=
8)
5000
HS
3.1802
1.1975
e
+
000
1.4046
e
+
000
1.2670
e
+
000
5.0691
e
-
002
IHS
3.1454
1.1667
e
+
000
1.3149
e
+
000
1.2283
e
+
000
3.8183
e
-
002
NGHS
3.1735
1.0260
e
+
000
1.1026
e
+
000
1.0445
e
+
000
1.8699
e
-
002
IGHS
4.7435
1.0200
e
+
000
1.0392
e
+
000
1.0277
e
+
000
5.3827
e
-
003Table4 gives the comparison of the results obtained by the IGHS, against the other three methods including the HS [15], the IHS [26], and the NGHS [20], and the best performance is reported in boldface. The terms “ACT” and “Std” stand for average computation time and standard deviation, respectively. According to the results, we can see that the IGHS performs better most of the time. To be precise, the values of Worst, Mean, and Std obtained by the IGHS are smaller than those obtained by the other three methods for all problems. Moreover, the IGHS can obtain five best objective function values except Example 1b (N
=
8). Therefore, the IGHS has exhibited stronger convergence and stability than the other three HSs. For Example 1b (N
=
8), the IHS, NGHS, and IGHS can yield the same best objective function value which is equal to 1.7119
e
-
003, and the remaining results obtained by the IGHS are better than the other three HSs.Additionally, Mann-WhitneyU test [27, 28], also known as “Mann-Whitney Wilcoxon test,” is used to ensure a statistical significant difference between the IGHS and any of the other three HSs. In other words, the Mann-Whitney U test acts as an efficient nonparametric rank-based test to identify a difference between populations. Moreover, the test statistic U is expressed as follows:(19)
U
1
=
R
1
-
n
1
n
1
+
1
2
,
(20)
U
2
=
R
2
-
n
2
n
2
+
1
2
,where n
1 and n
2 denote the sizes of sample 1 and sample 2, respectively. R
1 and R
2 denote the sums of the ranks in sample 1 and sample 2, respectively, U
1 represents the number of sample 1 observations beaten by sample 2 observations, and U
2 represents the number of sample 2 observations beaten by sample 1 observations. By combining (19) and (20), the sum of U
1 and U
2 is calculated as follows:(21)
U
1
+
U
2
=
R
1
-
n
1
n
1
+
1
2
+
R
2
-
n
2
n
2
+
1
2
.Based on the known conditionsR
1
+
R
2
=
n
(
n
+
1
)
/
2 and n
=
n
1
+
n
2, we can easily obtain the simplified form of (21) as follows:(22)
U
1
+
U
2
=
n
1
n
2
.To make the parametersR
1 and R
2 easy to understand, one simple and interesting example is provided as follows: a sample of 6 tortoises (sample 1) and 6 hares (sample 2) are collected and made run in a race. The order in which they reach the finishing post (their rank order, from first to last) is as follows: T H H H H H T T T T T H, where T represents a tortoise and H represents a hare.A direct way: we take each tortoise in turn and count the number of hares it is beaten by (lower rank), getting 0, 5, 5, 5, 5, and 5, which meansU
1
=
0
+
5
+
5
+
5
+
5
+
5
=
25. In the meantime, we could take each hare in turn and count the number of tortoises it is beaten by. In this situation, we get 1, 1, 1, 1, 1, and 6. So U
2
=
6
+
1
+
1
+
1
+
1
+
1
=
11. It is clear that the sum of these two values for U is 36, which is 6
×
6.An indirect way: the sum of the ranks achieved by the tortoises isR
1
=
1
+
7
+
8
+
9
+
10
+
11
=
46, and thus U
1
=
46
-
(
6
×
7
)
/
2
=
46
-
21
=
25 (as in (19)). The sum of the ranks achieved by the hares is R
2
=
2
+
3
+
4
+
5
+
6
+
12
=
32, and thus U
2
=
32
-
21
=
11 (as in (20)).In order to compare the IGHS with the other three HSs in a statistical way, three groups of Mann-WhitneyU tests are executed, and they are (
U
H
S
,
U
I
G
H
S
), (
U
I
H
S
,
U
I
G
H
S
), and (
U
N
G
H
S
,
U
I
G
H
S
), respectively. For each problem, 20 independent experiments are carried out; therefore, n
o
t
h
e
r
=
n
I
G
H
S
=
20 and U
o
t
h
e
r
+
U
I
G
H
S
=
n
o
t
h
e
r
×
n
I
G
H
S
=
400, where the subscript o
t
h
e
r denotes any of the other three HSs including the HS, IHS, and NGHS. The parameters of the four HSs are exactly the same as those of the aforementioned section, and Table 5 lists the results of Mann-Whitney U tests.Table 5
Mann-WhitneyU test result obtained using four methods.
Problem/algorithm
U
H
S
U
I
H
S
U
N
G
H
S
Example1a (N
=
5)
400
400
380
Example1a (N
=
8)
400
400
390
Example1b (N
=
5)
400
177
298
Example1b (N
=
8)
400
355
257
Example2a (N
=
8)
400
400
395
Example2b (N
=
8)
400
400
347From Table5, it is evident that the IGHS completely dominates the HS for solving all problems, because the values of U
H
S are all equal to 400. In other words, the IGHS has never been beaten by the HS. Furthermore, the values of U
N
G
H
S are significantly greater than 200, which indicates that the IGHS is also superior to the NGHS. In addition, the values of U
I
H
S are significantly greater than 200 for all problems except Example 1b (N
=
5). Thus, the IGHS performs better than the IHS on solving five of the six problems. Regarding Example 1b (N
=
5), the U
I
H
S value is equal to 177, which is smaller than but close to 200. Thus, the IHS is a little better than the IGHS for Example 1b (N
=
5).Figure8 shows the average convergence curves of four methods over 20 runs for six problems. It is clear that the HS has the poorest performance, and its convergence rate is the slowest. In addition, both the IHS and the NGHS converge faster than the HS but slower than the IGHS. Obviously, the IGHS has the fastest convergence rate in each case. For Example 1a (
N
=
5
) and (
N
=
8
), it can achieve lower levels than the other three HSs. With respect to Example 2b (
N
=
8
), the NGHS and the IGHS can converge to comparable levels, which are lower than those of the HS and the IHS. Moreover, the IHS, the NGHS, and the IGHS can reach comparable levels for the remaining cases. Strictly, the IGHS still outperforms the other three HSs for these cases according to Table 4.Figure 8
Average convergence curves of four methods on solving six problems.
(a)
Example1a (N
=
5
)
(b)
Example1a (N
=
8
)
(c)
Example1b (N
=
5
)
(d)
Example1b (N
=
8
)
(e)
Example2a (N
=
8
)
(f)
Example2b (N
=
8
)In Example1, the standard deviation of the Gaussian noise is so small (σ
=
0.001) as it is equivalent to considering a negligible noise, while in Example 2 the simulations are carried out without noise. To investigate the robustness of the proposed algorithm with respect to the additive noise, different standard deviations (σ) of the Gaussian noise are considered when solving a problem of system identification. On the other hand, signal to noise ratio (SNR) is a measure used in science and engineering that compares the level of a desired signal to the level of background noise. In this paper, the amplitudes of signal and noise are measured, and thus the SNR value can be calculated by S
N
R
d
B
=
20
l
o
g
10
(
A
s
i
g
n
a
l
/
A
n
o
i
s
e
), where A represents root mean square (RMS) amplitude. Given a standard deviation for a problem, a corresponding SNR value can be determined. Therefore, the SNR values related to various standard deviations are reported in Table 6.Table 6
The SNR values (in dB form) corresponding to various standard deviations.
Problem/algorithm
σ
=
0.001
σ
=
0.002
σ
=
0.005
σ
=
0.01
σ
=
0.02
σ
=
0.05
σ
=
0.1
σ
=
0.2
Example1a (N
=
5)
49.88
44.812
36.005
28.965
23.873
16.266
10.286
3.8217
Example1a (N
=
8)
49.567
43.463
35.865
29.496
24.641
15.688
10.307
3.911
Example1b (N
=
5)
51.415
46.475
37.818
31.1
25.635
17.753
12.19
5.296
Example1b (N
=
8)
51.448
44.949
37.822
31.875
25.634
19.216
11.73
5.2978
Example2a (N
=
8)
89.077
65.084
58.361
51.488
48.062
38.586
32.856
23.167
Example2b (N
=
8)
77.16
70.793
63.382
56.586
51.653
43.103
36.481
31.221From Table6, it is clear that the SNR value decreases as the standard deviation of the Gaussian noise increases. Therefore, a Gaussian noise with larger σ value will exert more serious effects on signal than the one with smaller σ value, which is harmful to the extraction of signal. In addition, a plot of MSE versus SNR is given in Figure 9.Figure 9
The mean square error (MSE) versus the signal to noise ratio (SNR).
(a)
Example1a (N
=
5
)
(b)
Example1a (N
=
8
)
(c)
Example1b (N
=
5
)
(d)
Example1b (N
=
8
)
(e)
Example2a (N
=
8
)
(f)
Example2b (N
=
8
)In Figures9(a), 9(b), 9(c), and 9(d), the fluctuations of MSE are minor in the range 20
≤
S
N
R
d
B
≤
55. In Figure 9(f), the fluctuations of MSE are minor in the range 40
≤
S
N
R
d
B
≤
80. In Figure 9(e), the fluctuations of MSE are slightly larger than those of the above five examples. Moreover, the MSE values tend to get larger when SNR decreases in most cases. In this experiment, the standard deviation σ varies in a wide range [
0.001,0.2
], the MSE values remain low levels for most σ values, and the large fluctuations of MSE do not happen until σ increases to very high levels (such as σ
=
0.1,0.2). Therefore, the robustness of IGHS with respect to the additive noise is acceptable to some extent.Various standard deviations of the Gaussian noise will have different effects on the convergence of HS, IHS, NGHS, and IGHS for six problems. To testify and compare the convergence of the proposed algorithm and the other three HSs, the average MSEs associated with the standard deviation are listed in Table7.Table 7
Effects of various standard deviations on the average MSEs of HS, IHS, NGHS, and IGHS.
Problem
NI
Algorithm
σ
=
0.002
σ
=
0.005
σ
=
0.01
σ
=
0.02
σ
=
0.05
σ
=
0.1
σ
=
0.2
Example1a (N
=
5)
2000
HS
0.031915
0.026339
0.021641
0.020161
0.031606
0.030344
0.073344
IHS
0.023133
0.01932
0.018847
0.017164
0.025531
0.02863
0.06602
NGHS
0.012373
0.01317
0.0080251
0.0087471
0.017701
0.015685
0.05541
IGHS
0.0090595
0.0081459
0.0063491
0.0060281
0.0082031
0.012456
0.05064
Example1a (N
=
8)
5000
HS
0.04305
0.053554
0.044272
0.061806
0.052752
0.056594
0.072702
IHS
0.040038
0.05833
0.039185
0.050786
0.059166
0.044153
0.065037
NGHS
0.01394
0.017639
0.011319
0.016349
0.017154
0.022171
0.034087
IGHS
0.0042096
0.0054617
0.0032124
0.004798
0.0065079
0.01026
0.021561
Example1b (N
=
5)
2000
HS
0.0020816
0.002023
0.0021274
0.0023366
0.0048039
0.011212
0.044087
IHS
0.0016932
0.0016984
0.001846
0.0020514
0.0042728
0.010803
0.043885
NGHS
0.0016883
0.0016907
0.0018335
0.0019744
0.004277
0.010804
0.043811
IGHS
0.0016881
0.0016894
0.0018321
0.0019734
0.0042728
0.010799
0.043809
Example1b (N
=
8)
5000
HS
0.0021598
0.002027
0.0018595
0.0023973
0.0035314
0.011308
0.042578
IHS
0.0017257
0.0017276
0.0017401
0.0022161
0.0034728
0.011198
0.042543
NGHS
0.0017311
0.0017232
0.0017334
0.0022162
0.0034424
0.011172
0.0425
IGHS
0.0017063
0.0017185
0.0017295
0.0022099
0.0034249
0.011163
0.042491
Example2a (N
=
8)
5000
HS
2.8627
3.9382
3.7033
3.2856
3.5386
5.1986
3.1728
IHS
2.8343
3.9123
3.6921
3.2648
3.5163
5.1678
3.1443
NGHS
2.6651
3.7442
3.5112
3.0952
3.3575
4.9886
2.9954
IGHS
2.6546
3.7404
3.4886
3.0747
3.3407
4.9715
2.9847
Example2b (N
=
8)
5000
HS
1.2711
1.2759
1.2776
1.2868
1.292
1.2358
1.3175
IHS
1.2339
1.236
1.2508
1.2473
1.242
1.2165
1.2836
NGHS
1.0437
1.0418
1.0432
1.0411
1.0456
1.0064
1.0917
IGHS
1.0276
1.0264
1.0294
1.029
1.0219
0.99213
1.0689The best results are marked in bold in Table7. For any problem with various standard deviations, IGHS can always obtain the smallest MSEs among four algorithms, indicating that IGHS has stronger convergence and stability than the other three HSs. Overall, IGHS has exhibited more desirable robustness to the Gaussian noise than the other three HSs for all six problems.


## 5.1. Example 1

The first example is the highly nonlinear discrete-time rational system [7], and its mathematical model is stated as follows:(15)
d
n
=
0.3
d
2
n
-
1
+
0.8
x
n
-
1
+
0.6
d
n
-
2
1
+
x
2
n
-
1
+
d
2
n
-
1
.Here, two types of input signals [7] are considered. The first is defined as Example 1a whose input is a random signal uniformly produced in [
-
1,1
]. The second is defined as Example 1b whose input is a cosine signal x
[
n
]
=
0.8
c
o
s
(
(
π
/
9
)
n
). For both Examples 1a and 1b, the measurement noise r
[
n
] is always supposed to be a Gaussian noise of N
(
0,0.001
).In this experiment, the IGHS is used for Example1a with N
=
5. Moreover, the IGHS parameters are set as follows: the maximal mutation rate p
¯
m
=
0.1, the minimal mutation rate p
_
m
=
0.1, harmony memory size H
M
S
=
5, the number of improvisations N
I
=
2000, and sampling number T
=
100. Moreover, T is actually the total number of output samples (actual system output or estimated output). According to (14), there are T sampling values for actual system output y
^
[
n
] (1
≤
n
≤
T). After implementing the IGHS, we can obtain the comparison between the estimated output y
[
n
] and the actual output y
^
[
n
] in Figure 2.Figure 2
Comparison of actual system output and Volterra model output for Example1a (N
=
5).From Figure2, it can be seen that the estimated output y
[
n
] approximates the actual output y
^
[
n
] well. In addition, the minimal mean square error (MSE) obtained by the IGHS is equal to 5.4594
e
-
003, which is a satisfactory result. Besides, we use Example 1a with N
=
8 to investigate the performance of the IGHS for second-order Volterra filter model with large memory size. The parameters of the IGHS are the same as those used for Example 1a (N
=
5) except N
I, and it is set to 5000 in this experiment. By performing the IGHS, we can obtain the comparison between the estimated output y
[
n
] and the actual output y
^
[
n
] in Figure 3.Figure 3
Comparison of actual system output and Volterra model output for Example1a (N
=
8).It is clear from Figure3 that a satisfactory approximation can be attained by utilizing the IGHS. Additionally, the minimal MSE yielded by the IGHS is equal to 3.7074
e
-
003 for Example 1a with N
=
8, which provides better modeling capacity.In addition to random signal, another testing input signalx
[
n
]
=
0.8
c
o
s
(
(
π
/
9
)
n
) is used to investigate the performance of second-order Volterra filter model using the IGHS. For Example 1b with N
=
5 and N
=
8, the IGHS parameters are the same as those for Example 1a, and Figures 4 and 5 display the comparisons of results for Example 1b with N
=
5 and N
=
8. As expected, the difference between the estimated output y
[
n
] and the actual output y
^
[
n
] is small in each case. Additionally, two satisfying MSEs can be obtained, respectively, for N
=
5 and N
=
8 by carrying out the IGHS, and they are equal to 1.71809876
e
-
003 and 1.7119
e
-
003, respectively.Figure 4
Comparison of actual system output and Volterra model output for Example1b (N
=
5).Figure 5
Comparison of actual system output and Volterra model output for Example1b (N
=
8).


## 5.2. Example 2

The second example is the real heat exchanger and its mathematical model [7, 25] is given by(16)
w
n
=
x
n
-
1.3228
x
2
n
+
0.7671
x
3
n
-
2.1755
x
4
n
,
(17)
d
n
=
-
6.5306
z
-
1
+
5.5652
z
-
2
1
-
1.608
z
-
1
+
0.6385
z
-
2
w
n
, where x
[
n
] is the process input denoting the flow rate and is constrained by the range of [
0,1
], w
[
n
] stands for the static nonlinearity, and d
[
n
] stands for the process output temperature. By combining (16) and (17), the simplified mathematical model can be expressed as the following difference equation:(18)
d
n
=
1.608
d
n
-
1
-
0.6385
d
n
-
2
-
6.5306
w
n
-
1
+
5.5652
w
n
-
2
.Here, two types of input signals [7, 25] are considered. The first is defined as Example 2a whose input is randomly generated in the range [
0.1,0.9
]. The second is defined as Example 2b whose input is a sine signal x
[
n
]
=
0.4
s
i
n
(
(
π
/
6
)
n
)
+
0.5. According to [7, 25], the measurement noise is excluded from the problem model; thus, we have r
[
n
]
=
0.In this experiment, the memory sizeN is set to 8 for second-order Volterra filter model. To fulfill the physical input requirement, the modeling input x
[
n
] is randomly generated in the range [
0.1,0.9
] and the testing input is then set to x
[
n
]
=
0.4
s
i
n
(
(
π
/
6
)
n
)
+
0.5. Moreover, the IGHS parameters used for Example 2 are the same as those of Example 1. By implementing the IGHS, the outputs of second-order Volterra model are shown in Figures 6 and 7, respectively, for the modeling input and testing input. Based on the careful observations on experimental results, we can confirm that the satisfactory approximation is attained in each case.Figure 6
Comparison of actual system output and Volterra model output for Example2a (N
=
8).Figure 7
Comparison of actual system output and Volterra model output for Example2b (N
=
8).


## 5.3. Comparison of the IGHS with the Other Three HSs

In this paper, four methods consisting of the HS [15], IHS [26], NGHS [20], and IGHS are used to solve the above problems with different input signals and system memory sizes. The parameters of these four methods are set as follows. For the HS, harmony memory considering rate H
M
C
R
=
0.95, pitch adjusting rate P
A
R
=
0.3, and bandwidth b
w
=
0.01. For the IHS, H
M
C
R
=
0.95, the minimal pitch adjusting rate P
A
R
m
i
n
=
0.35, the maximal pitch adjusting rate P
A
R
m
a
x
=
0.99, the minimal bandwidth b
w
m
i
n
=
1
0
-
5, and the maximal bandwidth b
w
m
a
x
=
0.1. For the NGHS, mutation rate p
m
=
5
×
1
0
-
2. For the IGHS, the minimal mutation rate p
_
m
=
0.01, the maximal mutation rate p
¯
m
=
0.1, and N
I
0
=
3
N
I
/
4. Note that all the methods are executed under identical conditions (harmony memory size and the number of improvisations). More specifically, the harmony memory size (H
M
S) of each method is set to 5, and the N
I (Number of improvisations) value of each method is set to 2000 and 5000, respectively, for the system memory sizes N
=
5 and N
=
8. With respect to second-order Volterra filter model, sampling number is set to T
=
100, and the range of each variable of the kernel vector is set as [
-
1,1
]. Matlab 7.0 is used to execute the above procedure under the environment of Intel(R) Core(TM) i5-2410M CPU @ 2.30 GHz. 20 independent runs are performed in each case, and the optimization results are presented in Table 4.Table 4
Comparison of HS, IHS, NGHS, and IGHS on six problems.
Problem
NI
Algorithm
ACT
Best
Worst
Mean
Std
Example1a (N
=
5)
2000
HS
1.659
1.1722
e
-
002
3.5945
e
-
002
1.9427
e
-
002
8.0506
e
-
003
IHS
1.7011
8.4989
e
-
003
2.6525
e
-
002
1.6498
e
-
002
5.9693
e
-
003
NGHS
1.781
5.9558
e
-
003
3.0649
e
-
002
1.2990
e
-
002
7.8612
e
-
003
IGHS
2.6365
5.4594
e
-
003
7.9210
e
-
003
5.9763
e
-
003
6.6092
e
-
004
Example1a (N
=
8)
5000
HS
4.3975
2.9721
e
-
002
9.3845
e
-
002
5.4468
e
-
002
1.5500
e
-
002
IHS
4.4121
3.1554
e
-
002
8.1672
e
-
002
5.1982
e
-
002
1.4302
e
-
002
NGHS
4.4334
5.9619
e
-
003
4.9489
e
-
002
1.5350
e
-
002
9.0064
e
-
003
IGHS
6.6589
3.7074
e
-
003
1.0237
e
-
002
5.6536
e
-
003
1.7608
e
-
003
Example1b (N
=
5)
2000
HS
0.57743
1.72096501
e
-
003
4.06941185
e
-
003
2.15386021
e
-
003
5.99250994
e
-
004
IHS
0.56219
1.71809880
e
-
003
1.82397928
e
-
003
1.72847844
e
-
003
2.71096253
e
-
005
NGHS
0.56886
1.71809919
e
-
003
1.83200101
e
-
003
1.72490455
e
-
003
2.54479841
e
-
005
IGHS
0.84918
1.71809876
e
-
003
1.71833108
e
-
003
1.71811926
e
-
003
5.18365165
e
-
008
Example1b (N
=
8)
5000
HS
2.6544
1.7202
e
-
003
2.3210
e
-
003
1.8309
e
-
003
1.5559
e
-
004
IHS
2.651
1.7119
e
-
003
1.7721
e
-
003
1.7194
e
-
003
1.8286
e
-
005
NGHS
2.6787
1.7119
e
-
003
1.8486
e
-
003
1.7188
e
-
003
3.0568
e
-
005
IGHS
4.0572
1.7119
e
-
003
1.7119
e
-
003
1.7119
e
-
003
1.4972
e
-
009
Example2a (N
=
8)
5000
HS
5.3252
3.8006
e
+
000
4.0102
e
+
000
3.8879
e
+
000
5.2382
e
-
002
IHS
5.2646
3.7901
e
+
000
3.9328
e
+
000
3.8548
e
+
000
4.3294
e
-
002
NGHS
5.3327
3.6604
e
+
000
3.7327
e
+
000
3.6760
e
+
000
1.6444
e
-
002
IGHS
7.9408
3.6560
e
+
000
3.6628
e
+
000
3.6581
e
+
000
1.9936
e
-
003
Example2b (N
=
8)
5000
HS
3.1802
1.1975
e
+
000
1.4046
e
+
000
1.2670
e
+
000
5.0691
e
-
002
IHS
3.1454
1.1667
e
+
000
1.3149
e
+
000
1.2283
e
+
000
3.8183
e
-
002
NGHS
3.1735
1.0260
e
+
000
1.1026
e
+
000
1.0445
e
+
000
1.8699
e
-
002
IGHS
4.7435
1.0200
e
+
000
1.0392
e
+
000
1.0277
e
+
000
5.3827
e
-
003Table4 gives the comparison of the results obtained by the IGHS, against the other three methods including the HS [15], the IHS [26], and the NGHS [20], and the best performance is reported in boldface. The terms “ACT” and “Std” stand for average computation time and standard deviation, respectively. According to the results, we can see that the IGHS performs better most of the time. To be precise, the values of Worst, Mean, and Std obtained by the IGHS are smaller than those obtained by the other three methods for all problems. Moreover, the IGHS can obtain five best objective function values except Example 1b (N
=
8). Therefore, the IGHS has exhibited stronger convergence and stability than the other three HSs. For Example 1b (N
=
8), the IHS, NGHS, and IGHS can yield the same best objective function value which is equal to 1.7119
e
-
003, and the remaining results obtained by the IGHS are better than the other three HSs.Additionally, Mann-WhitneyU test [27, 28], also known as “Mann-Whitney Wilcoxon test,” is used to ensure a statistical significant difference between the IGHS and any of the other three HSs. In other words, the Mann-Whitney U test acts as an efficient nonparametric rank-based test to identify a difference between populations. Moreover, the test statistic U is expressed as follows:(19)
U
1
=
R
1
-
n
1
n
1
+
1
2
,
(20)
U
2
=
R
2
-
n
2
n
2
+
1
2
,where n
1 and n
2 denote the sizes of sample 1 and sample 2, respectively. R
1 and R
2 denote the sums of the ranks in sample 1 and sample 2, respectively, U
1 represents the number of sample 1 observations beaten by sample 2 observations, and U
2 represents the number of sample 2 observations beaten by sample 1 observations. By combining (19) and (20), the sum of U
1 and U
2 is calculated as follows:(21)
U
1
+
U
2
=
R
1
-
n
1
n
1
+
1
2
+
R
2
-
n
2
n
2
+
1
2
.Based on the known conditionsR
1
+
R
2
=
n
(
n
+
1
)
/
2 and n
=
n
1
+
n
2, we can easily obtain the simplified form of (21) as follows:(22)
U
1
+
U
2
=
n
1
n
2
.To make the parametersR
1 and R
2 easy to understand, one simple and interesting example is provided as follows: a sample of 6 tortoises (sample 1) and 6 hares (sample 2) are collected and made run in a race. The order in which they reach the finishing post (their rank order, from first to last) is as follows: T H H H H H T T T T T H, where T represents a tortoise and H represents a hare.A direct way: we take each tortoise in turn and count the number of hares it is beaten by (lower rank), getting 0, 5, 5, 5, 5, and 5, which meansU
1
=
0
+
5
+
5
+
5
+
5
+
5
=
25. In the meantime, we could take each hare in turn and count the number of tortoises it is beaten by. In this situation, we get 1, 1, 1, 1, 1, and 6. So U
2
=
6
+
1
+
1
+
1
+
1
+
1
=
11. It is clear that the sum of these two values for U is 36, which is 6
×
6.An indirect way: the sum of the ranks achieved by the tortoises isR
1
=
1
+
7
+
8
+
9
+
10
+
11
=
46, and thus U
1
=
46
-
(
6
×
7
)
/
2
=
46
-
21
=
25 (as in (19)). The sum of the ranks achieved by the hares is R
2
=
2
+
3
+
4
+
5
+
6
+
12
=
32, and thus U
2
=
32
-
21
=
11 (as in (20)).In order to compare the IGHS with the other three HSs in a statistical way, three groups of Mann-WhitneyU tests are executed, and they are (
U
H
S
,
U
I
G
H
S
), (
U
I
H
S
,
U
I
G
H
S
), and (
U
N
G
H
S
,
U
I
G
H
S
), respectively. For each problem, 20 independent experiments are carried out; therefore, n
o
t
h
e
r
=
n
I
G
H
S
=
20 and U
o
t
h
e
r
+
U
I
G
H
S
=
n
o
t
h
e
r
×
n
I
G
H
S
=
400, where the subscript o
t
h
e
r denotes any of the other three HSs including the HS, IHS, and NGHS. The parameters of the four HSs are exactly the same as those of the aforementioned section, and Table 5 lists the results of Mann-Whitney U tests.Table 5
Mann-WhitneyU test result obtained using four methods.
Problem/algorithm
U
H
S
U
I
H
S
U
N
G
H
S
Example1a (N
=
5)
400
400
380
Example1a (N
=
8)
400
400
390
Example1b (N
=
5)
400
177
298
Example1b (N
=
8)
400
355
257
Example2a (N
=
8)
400
400
395
Example2b (N
=
8)
400
400
347From Table5, it is evident that the IGHS completely dominates the HS for solving all problems, because the values of U
H
S are all equal to 400. In other words, the IGHS has never been beaten by the HS. Furthermore, the values of U
N
G
H
S are significantly greater than 200, which indicates that the IGHS is also superior to the NGHS. In addition, the values of U
I
H
S are significantly greater than 200 for all problems except Example 1b (N
=
5). Thus, the IGHS performs better than the IHS on solving five of the six problems. Regarding Example 1b (N
=
5), the U
I
H
S value is equal to 177, which is smaller than but close to 200. Thus, the IHS is a little better than the IGHS for Example 1b (N
=
5).Figure8 shows the average convergence curves of four methods over 20 runs for six problems. It is clear that the HS has the poorest performance, and its convergence rate is the slowest. In addition, both the IHS and the NGHS converge faster than the HS but slower than the IGHS. Obviously, the IGHS has the fastest convergence rate in each case. For Example 1a (
N
=
5
) and (
N
=
8
), it can achieve lower levels than the other three HSs. With respect to Example 2b (
N
=
8
), the NGHS and the IGHS can converge to comparable levels, which are lower than those of the HS and the IHS. Moreover, the IHS, the NGHS, and the IGHS can reach comparable levels for the remaining cases. Strictly, the IGHS still outperforms the other three HSs for these cases according to Table 4.Figure 8
Average convergence curves of four methods on solving six problems.
(a)
Example1a (N
=
5
)
(b)
Example1a (N
=
8
)
(c)
Example1b (N
=
5
)
(d)
Example1b (N
=
8
)
(e)
Example2a (N
=
8
)
(f)
Example2b (N
=
8
)In Example1, the standard deviation of the Gaussian noise is so small (σ
=
0.001) as it is equivalent to considering a negligible noise, while in Example 2 the simulations are carried out without noise. To investigate the robustness of the proposed algorithm with respect to the additive noise, different standard deviations (σ) of the Gaussian noise are considered when solving a problem of system identification. On the other hand, signal to noise ratio (SNR) is a measure used in science and engineering that compares the level of a desired signal to the level of background noise. In this paper, the amplitudes of signal and noise are measured, and thus the SNR value can be calculated by S
N
R
d
B
=
20
l
o
g
10
(
A
s
i
g
n
a
l
/
A
n
o
i
s
e
), where A represents root mean square (RMS) amplitude. Given a standard deviation for a problem, a corresponding SNR value can be determined. Therefore, the SNR values related to various standard deviations are reported in Table 6.Table 6
The SNR values (in dB form) corresponding to various standard deviations.
Problem/algorithm
σ
=
0.001
σ
=
0.002
σ
=
0.005
σ
=
0.01
σ
=
0.02
σ
=
0.05
σ
=
0.1
σ
=
0.2
Example1a (N
=
5)
49.88
44.812
36.005
28.965
23.873
16.266
10.286
3.8217
Example1a (N
=
8)
49.567
43.463
35.865
29.496
24.641
15.688
10.307
3.911
Example1b (N
=
5)
51.415
46.475
37.818
31.1
25.635
17.753
12.19
5.296
Example1b (N
=
8)
51.448
44.949
37.822
31.875
25.634
19.216
11.73
5.2978
Example2a (N
=
8)
89.077
65.084
58.361
51.488
48.062
38.586
32.856
23.167
Example2b (N
=
8)
77.16
70.793
63.382
56.586
51.653
43.103
36.481
31.221From Table6, it is clear that the SNR value decreases as the standard deviation of the Gaussian noise increases. Therefore, a Gaussian noise with larger σ value will exert more serious effects on signal than the one with smaller σ value, which is harmful to the extraction of signal. In addition, a plot of MSE versus SNR is given in Figure 9.Figure 9
The mean square error (MSE) versus the signal to noise ratio (SNR).
(a)
Example1a (N
=
5
)
(b)
Example1a (N
=
8
)
(c)
Example1b (N
=
5
)
(d)
Example1b (N
=
8
)
(e)
Example2a (N
=
8
)
(f)
Example2b (N
=
8
)In Figures9(a), 9(b), 9(c), and 9(d), the fluctuations of MSE are minor in the range 20
≤
S
N
R
d
B
≤
55. In Figure 9(f), the fluctuations of MSE are minor in the range 40
≤
S
N
R
d
B
≤
80. In Figure 9(e), the fluctuations of MSE are slightly larger than those of the above five examples. Moreover, the MSE values tend to get larger when SNR decreases in most cases. In this experiment, the standard deviation σ varies in a wide range [
0.001,0.2
], the MSE values remain low levels for most σ values, and the large fluctuations of MSE do not happen until σ increases to very high levels (such as σ
=
0.1,0.2). Therefore, the robustness of IGHS with respect to the additive noise is acceptable to some extent.Various standard deviations of the Gaussian noise will have different effects on the convergence of HS, IHS, NGHS, and IGHS for six problems. To testify and compare the convergence of the proposed algorithm and the other three HSs, the average MSEs associated with the standard deviation are listed in Table7.Table 7
Effects of various standard deviations on the average MSEs of HS, IHS, NGHS, and IGHS.
Problem
NI
Algorithm
σ
=
0.002
σ
=
0.005
σ
=
0.01
σ
=
0.02
σ
=
0.05
σ
=
0.1
σ
=
0.2
Example1a (N
=
5)
2000
HS
0.031915
0.026339
0.021641
0.020161
0.031606
0.030344
0.073344
IHS
0.023133
0.01932
0.018847
0.017164
0.025531
0.02863
0.06602
NGHS
0.012373
0.01317
0.0080251
0.0087471
0.017701
0.015685
0.05541
IGHS
0.0090595
0.0081459
0.0063491
0.0060281
0.0082031
0.012456
0.05064
Example1a (N
=
8)
5000
HS
0.04305
0.053554
0.044272
0.061806
0.052752
0.056594
0.072702
IHS
0.040038
0.05833
0.039185
0.050786
0.059166
0.044153
0.065037
NGHS
0.01394
0.017639
0.011319
0.016349
0.017154
0.022171
0.034087
IGHS
0.0042096
0.0054617
0.0032124
0.004798
0.0065079
0.01026
0.021561
Example1b (N
=
5)
2000
HS
0.0020816
0.002023
0.0021274
0.0023366
0.0048039
0.011212
0.044087
IHS
0.0016932
0.0016984
0.001846
0.0020514
0.0042728
0.010803
0.043885
NGHS
0.0016883
0.0016907
0.0018335
0.0019744
0.004277
0.010804
0.043811
IGHS
0.0016881
0.0016894
0.0018321
0.0019734
0.0042728
0.010799
0.043809
Example1b (N
=
8)
5000
HS
0.0021598
0.002027
0.0018595
0.0023973
0.0035314
0.011308
0.042578
IHS
0.0017257
0.0017276
0.0017401
0.0022161
0.0034728
0.011198
0.042543
NGHS
0.0017311
0.0017232
0.0017334
0.0022162
0.0034424
0.011172
0.0425
IGHS
0.0017063
0.0017185
0.0017295
0.0022099
0.0034249
0.011163
0.042491
Example2a (N
=
8)
5000
HS
2.8627
3.9382
3.7033
3.2856
3.5386
5.1986
3.1728
IHS
2.8343
3.9123
3.6921
3.2648
3.5163
5.1678
3.1443
NGHS
2.6651
3.7442
3.5112
3.0952
3.3575
4.9886
2.9954
IGHS
2.6546
3.7404
3.4886
3.0747
3.3407
4.9715
2.9847
Example2b (N
=
8)
5000
HS
1.2711
1.2759
1.2776
1.2868
1.292
1.2358
1.3175
IHS
1.2339
1.236
1.2508
1.2473
1.242
1.2165
1.2836
NGHS
1.0437
1.0418
1.0432
1.0411
1.0456
1.0064
1.0917
IGHS
1.0276
1.0264
1.0294
1.029
1.0219
0.99213
1.0689The best results are marked in bold in Table7. For any problem with various standard deviations, IGHS can always obtain the smallest MSEs among four algorithms, indicating that IGHS has stronger convergence and stability than the other three HSs. Overall, IGHS has exhibited more desirable robustness to the Gaussian noise than the other three HSs for all six problems.


## 6. Conclusions

In the present paper, we propose an improved NGHS algorithm, called IGHS, for identifying nonlinear discrete-time systems based on second-order Volterra model. The prime objective of the IGHS is to attain the most appropriate kernels of Volterra model so that the estimated output can track and characterize the actual output as much as possible. Furthermore, we use two examples with different input signals and system memory sizes to test the performance of the IGHS. Experimental results suggest that the IGHS performs well and is superior to the other three methods in most cases according to four criteria (“Best,” “Worst,” “Mean,” and “Std”) associated with the mean square error (MSE). Thus, the IGHS is a strong candidate for the identification of nonlinear system based on second-order Volterra model.



---

*Source: 3102845-2016-02-09.xml*