# Identification of Nonlinear Dynamic Systems Using Hammerstein-Type Neural Network

Authors: Hongshan Yu; Jinzhu Peng; Yandong Tang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Hammerstein model has been popularly applied to identify the nonlinear systems. In this paper, a Hammerstein-type neural network (HTNN) is derived to formulate the well-known Hammerstein model. The HTNN consists of a nonlinear static gain in cascade with a linear dynamic part. First, the Lipschitz criterion for order determination is derived. Second, the backpropagation algorithm for updating the network weights is presented, and the stability analysis is also drawn. Finally, simulation results show that HTNN identification approach demonstrated identification performances.


---

## Body

## 1. Introduction

Identification of nonlinear dynamic systems has been one of the most interesting research areas in engineering. Similar as the well-known Wiener model, Hammerstein model is also well known and the most widely used for modeling of various processes [1, 2], which comprises of a static nonlinear block preceding a dynamic linear one [3]. This type of model was called block-oriented nonlinear model [4]. Different from black-box models, the block-oriented model was regarded as gray-box model which has clear physical interpretation, and its steady-state part describes the gain of the system [5]. The major drawback of these classic identification techniques is that they could have many constraints for nonlinear systems with time-varying parameters and uncertainties.It has been proven that neural networks (NN) can generally approximate any smooth nonlinear function well [2, 3]. During the last decade, they have been successfully and efficiently applied to identify and control the nonlinear systems due to high adaptation, self-organization, and fast parallel processing of real-time information. However, for general neural networks, the model structures have no relation with the physical nature of the process and the model parameters have no physical interpretation, similar to black-box models [5]. Therefore, many works were focused on integrating NN with Wiener or Hammerstein models. As a result, the model structures and the model parameters of NN have clear physical interpretation while the well-developed classic identification methodologies would have adaptive abilities.There are two ways to integrate NN with Wiener or Hammerstein models. One is to use NN to formulate the system static nonlinearities of Wiener or Hammerstein models. Chen et al. [6] used a simple linear model to represent the dynamic part and a neural network to represent the nonlinear static part. The dynamic linear part was replaced by Laguerre filters and the nonlinear static part was described as a neural network [2]. Tötterman and Toivonen [7] used support vector regression to identify nonlinear Wiener systems. The linear block was expanded in terms of Laguerre or Kautz filters, and the static nonlinear block was determined using support vector machine regression. Saha et al. [8] developed a Weiner-type nonlinear black-box model for capturing dynamics of open loop stable MIMO nonlinear systems, where quadratic polynomials and NN were used for constructing the nonlinear output map. The models were then used in nonlinear model predictive control [8]. Al-Duwaish et al. [9] used a hybrid model consisting of a linear autoregressive moving average (ARMA) model and a NN to, respectively, represent the dynamic linear block and the static nonlinear element of Wiener model. The other one is to formulate the Wiener or Hammerstein model entirely by a multilayer NN; that is, the dynamic linear block and static nonlinear block of Wiener or Hammerstein model are both represented by the NN. Janczak [4] designed a NN to formulate the Hammerstein model, which was composed of one hidden layer with nonlinear nodes and one linear output node. Wu et al. [10] proposed a Hammerstein neural network compensator to identify a dynamic process. Peng and Dubay [11] proposed a Wiener-type neural network to identify nonlinear dynamic processes. However, most of those investigations did not give the system order determination or stability analysis.In this paper, a multilayer neural network is used to formulate the traditional Hammerstein model, which is called Hammerstein-type neural network (HTNN). The HTNN formulates the Hammerstein model with a nonlinear static block in cascade with a linear dynamic block, the weights in which are corresponding with the parameters of Hammerstein model. Then, an identification methodology based on HTNN is presented and applied to nonlinear SISO systems. To determine the order and weights of HTNN, the Lipschitz criterion and backpropagation (BP) algorithm are derived, respectively. Furthermore, the stability analysis is drawn. Finally, the proposed identification method is tested on several nonlinear plants.In [2, 6–9], neural networks were used to formulate the system static nonlinearities of Wiener or Hammerstein models, while in our design, a multilayer neural network is used to formulate the Hammerstein model entirely; that is, the dynamic linear block and static nonlinear block of Hammerstein model were both represented by the neural network. In this way, the parameters of Hammerstein model can be obtained by training HTNN using an adequate training algorithm. Moreover, the model order determination and convergence stability are also drawn.The rest of this paper is organized as follows. In Section2, Hammerstein model is described. The design of the HTNN for identification is given in Section 3, the Lipschitz criterion is presented for order determination, and the learning algorithm is analyzed. The convergence analysis of HTNN is given in Section 4. The simulation results are given in Section 5. The conclusions are drawn in Section 6.


## 2. Hammerstein Model

Many industrial processes can be described by Wiener model or Hammerstein models [1]. As shown in Figure 1, a general Hammerstein model consists of a nonlinear static block and a linear dynamic block. The nonlinear static block is given by
(1)
x
(
t
)
=
f
(
u
(
t
)
)
,
where u
(
t
) is the input variable; f
(
·
) represents the nonlinear component of the Hammerstein model; x
(
t
) is a nonmeasured intermediate variable that does not necessarily have a physical meaning.Figure 1
The generally Hammerstein model.And the linear dynamic block can be described as(2)
y
(
t
)
=
G
(
q
-
1
)
x
(
t
)
=
A
(
q
-
1
)
B
(
q
-
1
)
x
(
t
)
with
(3)
A
(
q
-
1
)
=
a
1
q
-
1
+
a
2
q
-
2
+
⋯
+
a
n
a
q
-
n
a
B
(
q
-
1
)
=
1
+
b
1
q
-
1
+
b
2
q
-
2
+
⋯
+
b
n
b
q
-
n
b
,
where y
(
t
) is the output variable; q
-
1 is the unit delay operator; n
a and n
b are the orders of the linear dynamics and generally n
b
≤
n
a.There are many methods to identify the model parameters of Wiener model or Hammerstein model, such as the least squares method or its recursive version [12], maximum likelihood methods [13], correlation methods [14], linear optimization methods [15], and nonlinear optimization methods [2]. In this study, a multilayer neural network is designed to formulate the Hammerstein model. Therefore, the parameters of the Hammerstein model can be directly obtained through training the multilayer neural network by the BP training algorithm.


## 3. System Identification Using Hammerstein-Type Neural Network

A general neural network can be regarded as a black-box model [5]. In this section, a Hammerstein-type neural network (HTNN) is designed to formulate the Hammerstein model entirely. To determine the minimal model order, an order determination method based on Lipschitz quotients [16] is utilized. And for weights updating, the BP algorithm is presented.


### 3.1. Hammerstein-Type Neural Network

As shown in Figure2, a multilayer neural network called Hammerstein-type neural network is designed, which consists of a nonlinear static element and a single linear node with two tapped delay lines forming the model of the linear dynamic element.Figure 2
The structure of Hammerstein-type neural network.In Figure2, u
(
t
) is the input variable, y
^
(
t
) is the output variable, and x
^
(
t
) is a nonmeasured intermediate variable that does not necessarily have a physical meaning. Λ
i, i
=
1,2
,
…
,
p, means the order of input u, the parameters a
^
i, b
^
j, and c
^
k mean the weights of HTNN, and q
-
1 is the unit delay operator. Giving an input signal u
(
t
), initial condition of x
(
t
) and y
(
t
) (t
≤
0), and initial values of a
^
i, b
^
j, and c
^
k, we can obtain the HTNN output y
^
(
t
) (t
>
0).As shown in Figure2, a polynomial function is used in the nonlinear static block, the output of hidden layer x
^
(
t
) can be expressed as
(4)
x
^
(
t
)
=
f
(
u
(
t
)
)
=
∑
k
=
1
p
‍
c
^
k
u
k
(
t
)
,
and the output y
^
(
t
) can be expressed as
(5)
y
^
(
t
)
=
-
b
^
1
y
^
(
t
-
1
)
-
⋯
-
b
^
n
b
y
^
(
t
-
n
b
)
+
a
^
1
x
^
(
t
-
1
)
+
⋯
+
a
^
n
a
x
^
(
t
-
n
a
)
=
-
∑
j
=
1
n
b
‍
b
^
j
y
^
(
t
-
j
)
+
∑
i
=
1
n
a
‍
a
^
i
x
^
(
t
-
i
)
,
where c
^
k for k
=
1
,
…
,
p, b
^
j for j
=
1
,
…
,
n
b, and a
^
i for i
=
1
,
…
,
n
a are the weights of HTNN, which are associated with the parameters of Hammerstein model in (1) and (2). As a result, the parameters of the Hammerstein model can be expressed directly as the weights of the dynamic neural network. Then, the identified Hammerstein model can be obtained by training HTNN using an adequate training algorithm.


### 3.2. Model Order Determination

Usually, before training a neural network, it is important to determine how many neurons are in the network. In the HTNN, the structure of linear dynamic part is determined byn
a and n
b. As the most popular statistical model selection criterion, Akaike information criterion (AIC) [17] is a common method to determine the model order. However, the AIC for model order determination is usually subjected to a complex optimization process. The Lipschitz quotients criterion, which was proposed by He and Asada [16], can be utilized to determine the model order by analyzing the input-output data. Peng and Dubay [11] introduced the Lipschitz quotients criterion to determine the model order in a Wiener neural network.Consider a general nonlinear SISO dynamic system, which can be described as(6)
y
(
t
)
=
g
(
y
(
t
-
1
)
,
…
,
y
(
t
-
n
b
)
,
u
(
t
-
1
)
,
…
,
u
(
t
-
n
a
)
)
,
where y
(
t
) and u
(
t
) are the output and input variables of the nonlinear dynamic system, and n
b and n
a are the true orders of the output and input, respectively; g
(
·
) is a nonlinear function assumed to be continuous and smooth.Rewriting (6) in a compact form gives
(7)
y
=
g
(
χ
1
,
χ
2
,
…
,
χ
n
)
,
where n
=
n
b
+
n
a is the number of input variables; let χ
=
[
χ
1
,
χ
2
,
…
,
χ
n
]
T. To reconstruct the nonlinear function g
(
·
) from the input-output data pairs [
χ
(
i
)
,
y
(
i
)
]
i
=
1
N
set, where N
set is the number of data sets used for the model order determination. Define the Lipschitz quotient L
i
j as
(8)
L
i
j
=
|
y
(
i
)
-
y
(
j
)
|
|
χ
(
i
)
-
χ
(
j
)
|
,
i
≠
j
,
where |
χ
(
i
)
-
χ
(
j
)
| is the distance of two points in the input space and |
y
(
i
)
-
y
(
j
)
| is the difference between g
(
χ
(
i
)
) and g
(
χ
(
j
)
). For data points with a small distance |
χ
(
i
)
-
χ
(
j
)
| between them, the Lipschitz quotient L
i
j
(
n
) can be rewritten as
(9)
L
i
j
(
n
)
=
|
δ
y
|
(
δ
χ
1
)
2
+
(
δ
χ
2
)
2
+
⋯
+
(
δ
χ
n
)
2
,
where δ
y
=
y
(
i
)
-
y
(
j
) and δ
χ
r
=
χ
r
(
i
)
-
χ
r
(
j
) for r
=
1
,
…
,
n and the subscript n in L
i
j
(
n
) represents the number of input variables in (7). According to the investigation in He and Asada [16], the values of L
i
j
(
n
) can be used as indicator when one or more input variables are missing or in the case when one or more redundant input variables are included. For instance, if one of input variables, χ
n, is missing from the input set, the Lipschitz quotient L
i
j
(
n
-
1
) will be considerably larger than L
i
j
(
n
) or even unbounded. On the contrary, when an input variable χ
n
+
1 is included and the Lipschitz quotients L
i
j
(
n
) and L
i
j
(
n
+
1
) are calculated, if χ
n
+
1 is a redundant input variable, there will be a slight difference between L
i
j
(
n
) and L
i
j
(
n
+
1
) but not significant.To avoid the effect of measurement noise, the following index [16] is utilized to determine an appropriate order:
(10)
L
(
n
)
=
(
n
∏
s
=
1
m
‍
L
(
n
)
(
s
)
)
1
/
m
,
where L
(
n
)
(
s
) is the sth largest Lipschitz quotient among all L
i
j
(
n
) with the n input variables (
χ
1
,
…
,
χ
n
), and parameter m is a positive number, usually selected to be m
=
0.01
N
set
~
0.02
N
set. For testing purposes, the stop criterion can be defined as [18]
(11)
|
L
(
n
+
1
)
-
L
(
n
)
|
max
⁡
(
1
,
|
L
(
n
)
|
)
<
ϵ
,
where ϵ
>
0 is a prespecified threshold.It should be noted that the number of nodes in the nonlinear layer,p, is chosen manually. This is because, from empirical experience, for most cases, 3
≤
p
≤
6 can be chosen [11].


### 3.3. Learning Algorithm

During the neural network learning procedures, the weights are updated by the BP algorithm. An error function can be defined as(12)
J
(
w
,
t
)
=
1
2
(
y
(
t
)
-
y
^
(
t
)
)
2
=
1
2
e
(
t
)
2
,
where e
(
t
)
=
y
(
t
)
-
y
^
(
t
) is the identification error; y
(
t
) and y
^
(
t
) are the actual system output and the neural network output, respectively. Let w be the weights, which consist of the parameters c
^
k, b
^
j, and a
^
i. Applying the steepest descent method, the optimization target is characterized to minimize the error function equation (12) with respect to the weights of the network. Consider
(13)
∂
J
∂
w
=
-
e
(
t
)
∂
y
^
(
t
)
∂
w
,
where w represents the element w.The general update rule is expressed as(14)
w
(
t
+
1
)
=
w
(
t
)
+
Δ
w
(
t
)
=
w
(
t
)
-
η
∂
J
∂
w
=
w
(
t
)
+
η
e
(
k
)
∂
y
^
(
t
)
∂
w
,
where η is the training rate.Considering the element ofy
^
(
t
-
j
) also being function of the weights a
^
i and b
^
j, the partial derivative of the intermediate variables y
^
(
t
) to linear dynamic block weights a
^
i, b
^
j and intermediate variables x
^
(
t
-
i
) can be calculated as
(15)
∂
y
^
(
t
)
∂
b
^
j
=
-
y
^
(
t
-
j
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
b
^
j
,
j
=
1
,
…
,
n
b
,
(16)
∂
y
^
(
t
)
∂
a
^
i
=
x
(
t
-
i
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
a
^
i
,
i
=
1
,
…
,
n
a
,
(17)
∂
y
^
(
t
)
∂
x
^
(
t
-
i
)
=
a
^
i
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
x
^
(
t
-
i
)
,
i
=
1
,
…
,
n
a
.The partial derivative of intermediate variablesx
^
(
t
) to weights c
^
k can be calculated as
(18)
∂
x
^
(
t
)
∂
c
^
k
=
u
k
(
t
)
,
k
=
1
,
…
,
p
.According to (17) and (18), the partial derivative of the HTNN output y
^
(
t
) to weights c
^
k can be calculated as
(19)
∂
y
^
(
t
)
∂
c
^
k
=
∑
i
=
1
n
a
‍
(
∂
y
^
(
t
)
∂
x
^
(
t
-
i
)
·
∂
x
^
(
t
-
i
)
∂
c
^
k
)
=
∑
i
=
1
n
a
‍
a
^
i
u
k
(
t
-
i
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
c
^
k
,
k
=
1
,
…
,
p
.According to (14), the update law of c
^
k, a
^
i, and b
^
j can be calculated as
(20)
a
^
i
(
t
+
1
)
=
a
^
i
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
a
^
i
,
i
=
1
,
…
,
n
a
,
(21)
b
^
j
(
t
+
1
)
=
b
^
j
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
b
^
j
,
j
=
1
,
…
,
n
b
,
(22)
c
^
k
(
t
+
1
)
=
c
^
k
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
c
^
k
,
k
=
1
,
…
,
p
,
where the partial derivatives in (20), (21), and (22) are shown in (16), (15), and (19), respectively. From the above analysis, the structure of the system identification using the HTNN is shown in Figure 3.Figure 3
Model identification using HTNN.


## 3.1. Hammerstein-Type Neural Network

As shown in Figure2, a multilayer neural network called Hammerstein-type neural network is designed, which consists of a nonlinear static element and a single linear node with two tapped delay lines forming the model of the linear dynamic element.Figure 2
The structure of Hammerstein-type neural network.In Figure2, u
(
t
) is the input variable, y
^
(
t
) is the output variable, and x
^
(
t
) is a nonmeasured intermediate variable that does not necessarily have a physical meaning. Λ
i, i
=
1,2
,
…
,
p, means the order of input u, the parameters a
^
i, b
^
j, and c
^
k mean the weights of HTNN, and q
-
1 is the unit delay operator. Giving an input signal u
(
t
), initial condition of x
(
t
) and y
(
t
) (t
≤
0), and initial values of a
^
i, b
^
j, and c
^
k, we can obtain the HTNN output y
^
(
t
) (t
>
0).As shown in Figure2, a polynomial function is used in the nonlinear static block, the output of hidden layer x
^
(
t
) can be expressed as
(4)
x
^
(
t
)
=
f
(
u
(
t
)
)
=
∑
k
=
1
p
‍
c
^
k
u
k
(
t
)
,
and the output y
^
(
t
) can be expressed as
(5)
y
^
(
t
)
=
-
b
^
1
y
^
(
t
-
1
)
-
⋯
-
b
^
n
b
y
^
(
t
-
n
b
)
+
a
^
1
x
^
(
t
-
1
)
+
⋯
+
a
^
n
a
x
^
(
t
-
n
a
)
=
-
∑
j
=
1
n
b
‍
b
^
j
y
^
(
t
-
j
)
+
∑
i
=
1
n
a
‍
a
^
i
x
^
(
t
-
i
)
,
where c
^
k for k
=
1
,
…
,
p, b
^
j for j
=
1
,
…
,
n
b, and a
^
i for i
=
1
,
…
,
n
a are the weights of HTNN, which are associated with the parameters of Hammerstein model in (1) and (2). As a result, the parameters of the Hammerstein model can be expressed directly as the weights of the dynamic neural network. Then, the identified Hammerstein model can be obtained by training HTNN using an adequate training algorithm.


## 3.2. Model Order Determination

Usually, before training a neural network, it is important to determine how many neurons are in the network. In the HTNN, the structure of linear dynamic part is determined byn
a and n
b. As the most popular statistical model selection criterion, Akaike information criterion (AIC) [17] is a common method to determine the model order. However, the AIC for model order determination is usually subjected to a complex optimization process. The Lipschitz quotients criterion, which was proposed by He and Asada [16], can be utilized to determine the model order by analyzing the input-output data. Peng and Dubay [11] introduced the Lipschitz quotients criterion to determine the model order in a Wiener neural network.Consider a general nonlinear SISO dynamic system, which can be described as(6)
y
(
t
)
=
g
(
y
(
t
-
1
)
,
…
,
y
(
t
-
n
b
)
,
u
(
t
-
1
)
,
…
,
u
(
t
-
n
a
)
)
,
where y
(
t
) and u
(
t
) are the output and input variables of the nonlinear dynamic system, and n
b and n
a are the true orders of the output and input, respectively; g
(
·
) is a nonlinear function assumed to be continuous and smooth.Rewriting (6) in a compact form gives
(7)
y
=
g
(
χ
1
,
χ
2
,
…
,
χ
n
)
,
where n
=
n
b
+
n
a is the number of input variables; let χ
=
[
χ
1
,
χ
2
,
…
,
χ
n
]
T. To reconstruct the nonlinear function g
(
·
) from the input-output data pairs [
χ
(
i
)
,
y
(
i
)
]
i
=
1
N
set, where N
set is the number of data sets used for the model order determination. Define the Lipschitz quotient L
i
j as
(8)
L
i
j
=
|
y
(
i
)
-
y
(
j
)
|
|
χ
(
i
)
-
χ
(
j
)
|
,
i
≠
j
,
where |
χ
(
i
)
-
χ
(
j
)
| is the distance of two points in the input space and |
y
(
i
)
-
y
(
j
)
| is the difference between g
(
χ
(
i
)
) and g
(
χ
(
j
)
). For data points with a small distance |
χ
(
i
)
-
χ
(
j
)
| between them, the Lipschitz quotient L
i
j
(
n
) can be rewritten as
(9)
L
i
j
(
n
)
=
|
δ
y
|
(
δ
χ
1
)
2
+
(
δ
χ
2
)
2
+
⋯
+
(
δ
χ
n
)
2
,
where δ
y
=
y
(
i
)
-
y
(
j
) and δ
χ
r
=
χ
r
(
i
)
-
χ
r
(
j
) for r
=
1
,
…
,
n and the subscript n in L
i
j
(
n
) represents the number of input variables in (7). According to the investigation in He and Asada [16], the values of L
i
j
(
n
) can be used as indicator when one or more input variables are missing or in the case when one or more redundant input variables are included. For instance, if one of input variables, χ
n, is missing from the input set, the Lipschitz quotient L
i
j
(
n
-
1
) will be considerably larger than L
i
j
(
n
) or even unbounded. On the contrary, when an input variable χ
n
+
1 is included and the Lipschitz quotients L
i
j
(
n
) and L
i
j
(
n
+
1
) are calculated, if χ
n
+
1 is a redundant input variable, there will be a slight difference between L
i
j
(
n
) and L
i
j
(
n
+
1
) but not significant.To avoid the effect of measurement noise, the following index [16] is utilized to determine an appropriate order:
(10)
L
(
n
)
=
(
n
∏
s
=
1
m
‍
L
(
n
)
(
s
)
)
1
/
m
,
where L
(
n
)
(
s
) is the sth largest Lipschitz quotient among all L
i
j
(
n
) with the n input variables (
χ
1
,
…
,
χ
n
), and parameter m is a positive number, usually selected to be m
=
0.01
N
set
~
0.02
N
set. For testing purposes, the stop criterion can be defined as [18]
(11)
|
L
(
n
+
1
)
-
L
(
n
)
|
max
⁡
(
1
,
|
L
(
n
)
|
)
<
ϵ
,
where ϵ
>
0 is a prespecified threshold.It should be noted that the number of nodes in the nonlinear layer,p, is chosen manually. This is because, from empirical experience, for most cases, 3
≤
p
≤
6 can be chosen [11].


## 3.3. Learning Algorithm

During the neural network learning procedures, the weights are updated by the BP algorithm. An error function can be defined as(12)
J
(
w
,
t
)
=
1
2
(
y
(
t
)
-
y
^
(
t
)
)
2
=
1
2
e
(
t
)
2
,
where e
(
t
)
=
y
(
t
)
-
y
^
(
t
) is the identification error; y
(
t
) and y
^
(
t
) are the actual system output and the neural network output, respectively. Let w be the weights, which consist of the parameters c
^
k, b
^
j, and a
^
i. Applying the steepest descent method, the optimization target is characterized to minimize the error function equation (12) with respect to the weights of the network. Consider
(13)
∂
J
∂
w
=
-
e
(
t
)
∂
y
^
(
t
)
∂
w
,
where w represents the element w.The general update rule is expressed as(14)
w
(
t
+
1
)
=
w
(
t
)
+
Δ
w
(
t
)
=
w
(
t
)
-
η
∂
J
∂
w
=
w
(
t
)
+
η
e
(
k
)
∂
y
^
(
t
)
∂
w
,
where η is the training rate.Considering the element ofy
^
(
t
-
j
) also being function of the weights a
^
i and b
^
j, the partial derivative of the intermediate variables y
^
(
t
) to linear dynamic block weights a
^
i, b
^
j and intermediate variables x
^
(
t
-
i
) can be calculated as
(15)
∂
y
^
(
t
)
∂
b
^
j
=
-
y
^
(
t
-
j
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
b
^
j
,
j
=
1
,
…
,
n
b
,
(16)
∂
y
^
(
t
)
∂
a
^
i
=
x
(
t
-
i
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
a
^
i
,
i
=
1
,
…
,
n
a
,
(17)
∂
y
^
(
t
)
∂
x
^
(
t
-
i
)
=
a
^
i
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
x
^
(
t
-
i
)
,
i
=
1
,
…
,
n
a
.The partial derivative of intermediate variablesx
^
(
t
) to weights c
^
k can be calculated as
(18)
∂
x
^
(
t
)
∂
c
^
k
=
u
k
(
t
)
,
k
=
1
,
…
,
p
.According to (17) and (18), the partial derivative of the HTNN output y
^
(
t
) to weights c
^
k can be calculated as
(19)
∂
y
^
(
t
)
∂
c
^
k
=
∑
i
=
1
n
a
‍
(
∂
y
^
(
t
)
∂
x
^
(
t
-
i
)
·
∂
x
^
(
t
-
i
)
∂
c
^
k
)
=
∑
i
=
1
n
a
‍
a
^
i
u
k
(
t
-
i
)
-
∑
m
=
1
n
b
‍
b
^
m
∂
y
^
(
t
-
m
)
∂
c
^
k
,
k
=
1
,
…
,
p
.According to (14), the update law of c
^
k, a
^
i, and b
^
j can be calculated as
(20)
a
^
i
(
t
+
1
)
=
a
^
i
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
a
^
i
,
i
=
1
,
…
,
n
a
,
(21)
b
^
j
(
t
+
1
)
=
b
^
j
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
b
^
j
,
j
=
1
,
…
,
n
b
,
(22)
c
^
k
(
t
+
1
)
=
c
^
k
(
t
)
+
η
·
e
^
(
t
)
∂
y
^
(
t
)
∂
c
^
k
,
k
=
1
,
…
,
p
,
where the partial derivatives in (20), (21), and (22) are shown in (16), (15), and (19), respectively. From the above analysis, the structure of the system identification using the HTNN is shown in Figure 3.Figure 3
Model identification using HTNN.


## 4. Convergence Analysis

In the training procedures of HTNN, a proper choice of training rateη is usually required in the update rules of (20)–(22). Too small η guarantees convergence but with slow training speed, while too big η leads to being unstable. In this section, the approach on selecting properly η is developed.A discrete type Lyapunov function can be defined as [19]
(23)
V
(
t
)
=
1
2
e
2
(
t
)
.Therefore, the time change of the Lyapunov function can be obtained as(24)
Δ
V
(
t
)
=
V
(
t
+
1
)
-
V
(
t
)
=
1
2
[
e
2
(
t
+
1
)
-
e
2
(
t
)
]
=
Δ
e
(
t
)
[
e
(
t
)
+
Δ
e
(
t
)
2
]
.The error difference can be represented as(25)
e
(
t
+
1
)
=
e
(
t
)
+
Δ
e
(
t
)
=
e
(
t
)
+
(
∂
e
(
t
)
∂
w
)
T
·
Δ
w
.From the update rule of (13), we have
(26)
Δ
w
=
-
η
·
e
(
t
)
·
∂
e
(
t
)
∂
w
=
η
·
e
(
t
)
·
∂
y
^
(
t
)
∂
w
.Theorem 1.
Assume thatη is the learning rate for the weights of HTNN and β
max
⁡ is defined as β
max
⁡
=
max
⁡
t
∥
β
(
t
)
∥, where β
(
t
)
=
∂
y
^
(
t
)
/
∂
w and ∥
·
∥ is the Euclidean norm. Then, the convergence is guaranteed if η is satisfied as
(27)
0
<
η
<
2
β
max
⁡
2
.Proof.
According to (23)-(24), Δ
V
(
t
) can be represented as
(28)
Δ
V
(
t
)
=
Δ
e
(
t
)
[
e
(
t
)
+
Δ
e
(
t
)
2
]
=
(
∂
e
(
t
)
∂
w
)
T
η
e
(
t
)
∂
y
^
(
t
)
∂
w
·
[
e
(
t
)
+
1
2
(
∂
e
(
t
)
∂
w
)
T
η
e
(
t
)
∂
y
^
(
t
)
∂
w
]
.
Substituting the equation ∂
e
(
t
)
/
∂
w
=
-
∂
y
^
(
t
)
/
∂
w into (28) yields
(29)
Δ
V
(
t
)
=
-
η
e
2
(
t
)
∥
∂
y
^
(
k
)
∂
w
∥
2
+
1
2
η
2
e
2
(
t
)
∥
∂
y
^
(
t
)
∂
w
∥
4
.
Define β
max
⁡
=
max
⁡
t
∥
β
(
t
)
∥, where β
(
t
)
=
∂
y
^
(
t
)
/
∂
w; we obtain
(30)
Δ
V
(
t
)
=
-
1
2
∥
β
(
t
)
∥
2
η
(
2
-
η
∥
β
(
t
)
∥
2
)
e
2
(
t
)
≤
-
1
2
∥
β
(
t
)
∥
2
η
(
2
-
η
β
max
⁡
2
)
e
2
(
t
)
.
Since V
(
t
)
>
0 for all time t, the convergence of training algorithm means Δ
V
(
t
)
<
0. Referring to (30), it implies that (27) is satisfied. It should be pointed that optimal convergence via maximum learning rate corresponds to η
*
=
1
/
β
max
⁡
2, which is the upper half of the limit in (27).


## 5. Simulation Examples

In this section, three examples are utilized to illustrate the HTNN identifier. The first is a Hammerstein process which is utilized to demonstrate the HTNN of fitting the Hammerstein model. The second and the third are a nonlinear dynamic system and an industry process, which are utilized to demonstrate the HTNN can be used to identify nonlinear dynamic systems and processes, respectively.The implementation procedure of the HTNN for nonlinear system identification is itemized as follows.Step 1.
Select the input-output variablesu
(
t
) and y
(
t
) and the structures (n
a and n
b) of the HTNN according to Section 3.2.Step 2.
Select the orderp of polynomial function; note that the value range of integer p is 3
≤
p
≤
6.Step 3.
Model the nonlinear system using HTNN as in (4) and (5), and train the HTNN using the training datasets according to (20)–(22) to obtain the weights a
^
i, b
^
j, and c
^
k.Step 4.
Calculate the error functionJ in (12). If J is less than a limiting value within a given number of training iterations, it means that the model can be accepted; otherwise go to Step 2 to reselect the value of integer p.Example 1.
Nonlinear Hammerstein model identification is as follows.
The following Hammerstein model is described as(31)
x
(
t
)
=
u
(
t
)
-
0.4
u
(
t
)
2
+
0.1
u
(
t
)
3
;
y
(
t
)
=
1.4138
y
(
t
-
1
)
-
0.6065
y
(
t
-
2
)
+
0.1044
x
(
t
-
1
)
+
0.0883
x
(
t
-
2
)
.
A total of 1000 data pairs of ani.i.d. uniform sequence within the limits [
-
0.5,0.5
] were generated to train the HTNN. In the order determination procedure, firstly, it gives the input χ
1
=
y
(
t
-
1
) and computes the Lipschitz quotient L
(
1,0
)
=
+
∞ and then sets χ
2
=
u
(
t
-
1
); the Lipschitz quotient L
(
1,1
)
=
85.55. For χ
3
=
y
(
t
-
2
) and χ
4
=
u
(
t
-
2
), the corresponding Lipschitz quotients L
(
2,1
)
=
6.312 and L
(
2,2
)
=
1.853 are decreased significantly. Then it sets χ
5
=
y
(
t
-
3
), χ
6
=
u
(
t
-
3
), and χ
7
=
y
(
t
-
4
), the corresponding Lipschitz quotients L
(
3,2
)
=
1.452, L
(
3,3
)
=
1.002, and L
(
4,3
)
=
0.9799 are not significantly different from L
(
2,2
), and the stop criterion equation (11) is satisfied, where the threshold is chosen as ϵ
=
0.5. Figure 4 shows the values of the Lipschitz quotients for different order; the best order of the system is (2,2); that is, n
a
=
2 and n
b
=
2. From (31), the true order of the system is (2,2).
In this case, the number of neurons in nonlinear static block is chosen asp
=
3. According to the convergence analysis, β
max
⁡
=
2.534; therefore, the learning rate can be chosen as η
=
0.16. The initial parameters in (21)–(23) are y
(
t
)
=
0, x
(
t
)
=
0, ∂
y
(
t
)
/
∂
a
^
i
=
0, ∂
y
(
t
)
/
∂
b
^
j
=
0, and ∂
y
(
t
)
/
∂
c
^
k
=
0 as t
≤
0. Then, the testing input signal u
(
t
)
=
0.5
sin
(
t
) is used to verify the identification performance of the HTNN. Figure 5 illustrates the output of the plant and the HTNN. And the mean square error (MSE) is 1.497
×
1
0
-
3 with 7 tunable parameters using the proposed HTNN.Figure 4
Values of the order determination.Figure 5
Outputs of the plant (solid) and HTNN (dotted).Example 2.
Nonlinear dynamic system identification is as follows.
The following process is a nonlinear dynamic process formulated in a discrete form as [20]
(32)
y
(
t
)
=
0.5
(
y
(
t
-
1
)
+
u
(
t
-
1
)
e
-
3
|
y
(
t
-
1
)
|
)
.
A total of 500 data pairs which are generated by inputs with a sinusoidal function1.05
sin
(
π
t
/
45
) are utilized to train the network.
Similar as Example1, the values of the Lipschitz quotients are L
(
1,0
)
=
+
∞, and L
(
1,1
)
=
41.45, L
(
2,1
)
=
2.761, L
(
3,1
)
=
1.596, L
(
2,2
)
=
1.389, L
(
3,2
)
=
1.157, L
(
3,3
)
=
0.9093, L
(
4,3
)
=
0.6388. The stop criterion ends at L
(
3,1
); it implies that the best order of the system is (3,1). Figure 6 shows the values of the Lipschitz quotients for different orders.
In this case, the number of neurons in linear dynamic block is chosen asn
a
=
3, n
b
=
1 and the number of neurons in nonlinear static block is chosen as p
=
3. According to the convergence analysis, β
max
⁡
=
2.108; therefore, the learning rate can be chosen as η
=
0.22. The testing input signal u
(
t
)
=
sin
(
2
π
t
)
/
10 is used to verify the identification performance of the HTNN.Figure 6
Values of the order determination.The proposed HTNN was compared to several neural network based identification methods, the Controllable-Canonical-Form-Based Recurrent Neurofuzzy Network (CReNN) [21] and the Dynamic Fuzzy Neural Network (DFNN) [22]. As for the DFNN model, it is a fuzzy model of three rules, where the two-dimensional input space was partitioned to three clusters and a Gaussian membership function was assigned to each cluster. Also, using least square (LS) algorithm, a standard Hammerstein model (SHM) with n
a
=
3, n
b
=
1, and p
=
3, which is shown in (1) and (2), was applied to identify the nonlinear system for comparison. As shown in Table 1, the results illustrate that the proposed HTNN has the least number of parameters with lower MSE value. Figure 7 shows the output of the plant, the standard Hammerstein model with least square algorithm (SHM-LS), and the HTNN identifier. It can be seen that the standard Hammerstein model with least square algorithm (SHM-LS) can not identify the complex nonlinear system well.Table 1
Comparison among several neural networks for identification.
Network
Parameters
MSE
HTNN
7
1.828
×
1
0
-
3
CReNN
51
4.784
×
1
0
-
2
DFNN
39
1.325
×
1
0
-
3
SHM-LS
7
0.2625Figure 7
Outputs of the plant (solid), HTNN (dotted), and SHM-LS (dash).Example 3.
Identification for an industry process is as follows.
The proposed HTNN was also used to identify a typical industry process, a continuous stirred tank reactor (CSTR). An irreversible first-order reaction is part of the CSTR, which has the dimensionless mass and energy balances. The system is described by the following equation [23, 24]:
(33)
d
C
A
d
t
=
Q
V
(
C
A
f
-
C
A
)
-
k
0
C
A
e
-
E
/
RT
d
T
d
t
=
Q
V
(
T
f
-
T
)
-
Δ
H
ρ
C
p
k
0
C
A
e
-
E
/
RT
+
U
A
V
ρ
C
p
(
T
c
-
T
)
,
where the description and value of each variable are given in Table 2. Two state variables of the model are the reactant concentration C
A and the reactor temperature T. The control objective is to control the reactant concentration C
A, through the manipulation of the coolant temperature T
c. Note that the reactor temperature T is not controlled for this simulation. Therefore, the output variable and manipulated variable are given by y
(
t
)
=
C
A and u
(
t
)
=
T
c, respectively. For practicability, the coolant temperature T
c is constrained to the ranges [
273.0,373.0
].Table 2
Model parameters of CSTR process.
Variables
Description
Values
Q
Volumetric flow rate
10 (L·min−1)
V
Reactor volume
10 (L)
ρ
Density of reaction mixture
1000 (g·L−1)
C
A
f
Feed concentration
1.0 (mol·L−1)
T
f
Feed temperature
350 (K)
C
p
Specific heat capacity
1.0 (J·g−1K−1)
Δ
H
Heat of reaction
-
1.0
×
1
0
5 (J·mol−1)
k
0
Arrhenius preexponential constant
5.33685
×
1
0
7 (min−1)
E
/
R
Activation energy/gas law constant
6000 (K)
U
A
Heat transfer term
5000 (J·min−1K−1)The above model is used to generate a series of input-output time-series data. The sampling time of the process measurements is set to 0.1 min. From Figure8, the stop criterion ends at L
(
2,1
)
=
2.333, indicating that the best order of the system is (2,1); that is, n
a
=
2 and n
b
=
1. The number of neurons in nonlinear static block p
=
4. According to the convergence analysis, β
max
⁡
=
2.856; therefore, the learning rate can be chosen as η
=
0.15, and the initial parameters in (21)–(23) are the same as Example 1. The nonlinear part of HTNN is sensitive to the data between -
1 and 1; however, the input values (T
c) are within [
273.0,373.0
]. Therefore, it is necessary to normalize the input signals as follows:
(34)
u
(
t
)
=
T
c
(
t
)
-
273.0
373.0
-
273.0
-
373.0
-
T
c
(
t
)
373.0
-
273.0
.Figure 8
Values of the order determination.The normalized data is then used to train the HTNN; the training results are shown in Figure9. The models were then tested on a set of data produced from an input with random amplitude; a common multilayer neural network with 5 hidden neurons [25] was also developed for comparison. As shown in Figure 9, the HTNN gives a good fit to this data.Figure 9
Training results, outputs of the plant (solid curve) and HTNN (dotted curve) and multilayer neural network (dash curve).


## 6. Conclusions

In this paper, by formulated Hammerstein model, a Hammerstein-type neural network was developed for identifying nonlinear SISO systems, where the weights are corresponding with the parameters of Hammerstein model. To determine the model order and the parameters of the Hammerstein model, the Lipschitz quotients and the backpropagation training algorithm were used to determine the model order and update the weights in the network, respectively, and the stability was also analyzed. The HTNN was tested on several nonlinear systems to demonstrate the identification performances.



---

*Source: 959507-2014-10-19.xml*