# Filtering Based Recursive Least Squares Algorithm for Multi-Input Multioutput Hammerstein Models

Authors: Ziyun Wang; Yan Wang; Zhicheng Ji
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

This paper considers the parameter estimation problem for Hammerstein multi-input multioutput finite impulse response (FIR-MA) systems. Filtered by the noise transfer function, the FIR-MA model is transformed into a controlled autoregressive model. The key-term variable separation principle is used to derive a data filtering based recursive least squares algorithm. The numerical examples confirm that the proposed algorithm can estimate parameters more accurately and has a higher computational efficiency compared with the recursive least squares algorithm.


---

## Body

## 1. Introduction

Parameter estimation is an important approach to model dynamical systems and has been widely used in estimating the parameters for nonlinear systems [1–3], deriving system identification methods [4–7], identifying state-space models [8, 9], and developing solutions for matrix equations [10–13]. For example, Dehghan and Hajarian discussed several solution methods for different matrix equations [14–16]. In the area of system control and modeling, Shi and Fang developed a Kalman filter based identification for systems with randomly missing measurements [17], gave output feedback stabilization [18], and presented a robust mixed H
2
/
H
∞ control of networked control systems [19].The least squares algorithm is a fundamental method [20–22] and many methods such as the iterative algorithm [23, 24] and the gradient algorithm [25] are widely used in the parameter estimation. In the field of Hammerstein system identification, several methods have been developed [26]. For example, a least squares based iterative algorithm and an auxiliary model based recursive least squares algorithm have been presented, respectively, for Hammerstein nonlinear ARMAX systems and Hammerstein output error systems [27, 28]; a Newton recursive algorithm and a Newton iterative algorithm for Hammerstein controlled autoregressive systems are presented in [29].Consider a multi-input multioutput (MIMO) Hammerstein finite impulse response (FIR) system depicted by(1)
y
(
t
)
=
B
(
z
)
u
-
(
t
)
+
N
(
z
)
v
(
t
)
,
where u
(
t
)
∶
=
[
u
1
(
t
)
,
u
2
(
t
)
,
…
,
u
r
(
t
)
]
T
∈
R
r is the nonlinear system input vector with zero mean and unit variances, y
(
t
)
∈
R
m is the measurement of x
(
t
)
∶
=
B
(
z
)
u
-
(
t
) but is corrupted by w
(
t
)
∶
=
N
(
z
)
v
(
t
), v
(
t
)
∈
R
m is the white noise vector with zero mean, and B
(
z
) and N
(
z
) are polynomials in the unit backward shift operator z
-
1
[
z
-
1
y
(
t
)
=
y
(
t
-
1
)
]:
(2)
B
(
z
)
∶
=
I
m
×
r
+
∑
i
=
1
n
b
‍
B
i
z
-
i
,
B
i
∈
R
m
×
r
,
N
(
z
)
∶
=
1
+
∑
j
=
1
n
d
‍
d
j
z
-
j
,
d
j
∈
R
1
.
It is obvious that the relation between sizes m and r would influence the model identification of this multi-input multioutput Hammerstein system. For example, the dimension of the output vector is not less than that of the input vector if m
⩾
r; otherwise, when m
<
r, the output size is smaller compared with that of the input vector. In this paper, we discuss the identification problem with m
⩾
r.The nonlinear block in the Hammerstein model is a linear combination of the known basisf
∶
=
(
f
1
,
f
2
,
…
,
f
r
):
(3)
u
-
(
t
)
=
f
(
u
(
t
)
)
∶
=
[
f
1
(
u
1
(
t
)
)
,
f
2
(
u
2
(
t
)
)
,
…
,
f
r
(
u
r
(
t
)
)
]
T
∈
R
r
,
where the superscript T denotes the matrix transpose. The function f
i
(
u
i
(
t
)
) in (3) is a nonlinear function of a known basis (
γ
1
,
γ
2
,
…
,
γ
n
c
):
(4)
f
i
(
u
i
(
t
)
)
∶
=
∑
j
=
1
n
c
‍
c
j
γ
j
(
u
i
(
t
)
)
,
where the coefficients c
∶
=
(
c
1
,
c
2
,
…
,
c
n
c
) are unknown. Substituting (4) into (3) yields
(5)
u
-
(
t
)
=
∑
i
=
1
n
c
‍
c
i
γ
i
(
u
(
t
)
)
,
γ
i
(
u
(
t
)
)
∶
=
[
γ
i
(
u
1
(
t
)
)
,
γ
i
(
u
2
(
t
)
)
,
…
,
γ
i
(
u
r
(
t
)
)
]
T
∈
R
r
.
Assume that u
(
t
)
=
0, y
(
t
)
=
0, and v
(
t
)
=
0 for t
⩽
0, and the orders n
b, n
c, and n
d are known but can be obtained by trial and error. In general, the orders of the Hammerstein model should be large when the nonlinear system is used in prediction; otherwise, the orders should be small if the system is applied for control.The objective of this paper is to estimate the unknown parameter matrices:B
i, c
i, d
i from the available input-output data {
u
(
t
)
,
y
(
t
)
} of the multivariable Hammerstein finite impulse response moving average (FIR-MA) models [30].Recently, the filtering idea has received much attention [31–33]. Xiao and Yue studied input nonlinear dynamical adjustment models and presented a recursive generalized least squares algorithm and a filtering based least squares algorithm by replacing the unknown terms in the information vectors with their estimates [34]. The overparameterization method in [34] leads to a redundant estimated product of the nonlinear systems and requires extra computation. Differing from the work in [30, 34, 35], this paper discusses the estimation problem of the MIMO Hammerstein systems using the data filtering idea and transfers the FIR-MA system to controlled autoregressive model by means of the key-term variable separate principle in [36–38]. The proposed algorithm used in this paper can extend to study parameter estimation problems of dual-rate/multirate sampled systems [39–42] and other linear or nonlinear systems [43–46].Briefly, the rest of this paper is recognized as follows. Section2 discusses a recursive least squares algorithm for the Hammerstein systems. Section 3 presents a filtering based recursive least squares algorithm by transferring an FIR-MA model to a controlled autoregressive model. Section 4 provides an illustrative example. Finally, some concluding remarks are offered in Section 5.


## 2. The MRLS Algorithm

For comparison, the MRLS algorithm is listed in Section2. Here we introduce some notations. t represents the current time in this paper and “A
=
:
X” or “X
∶
=
A” means that “A is defined as X”. The symbol I
m
×
r represents an identity matrix of size r followed by a null matrix of the last m
-
r rows when m
⩾
r and vice versa. The norm of a matrix (or a column vector) X is defined by ∥
X
∥
2
∶
=
tr
⁡
[
X
X
T
]; ⊗ denotes the Kronecker product or direct product: if A
=
[
a
i
j
]
∈
R
m
×
n, B
=
[
b
i
j
]
∈
R
p
×
q, then A
⊗
B
∶
=
[
a
i
j
B
]
∈
R
m
p
×
n
q; col[
X
] is supposed to be the vector formed by the column of the matrix X: if X
=
[
x
1
,
x
2
,
…
,
x
n
]
∈
R
m
×
n, then col[
X
]
∶
=
[
x
1
T
,
x
2
T
,
…
,
x
n
T
]
T
∈
R
m
n.From (1)–(5), the intermediate variables x
(
t
) and w
(
t
) and output of the system y
(
t
) can be expressed as
(6)
x
(
t
)
=
[
I
m
×
r
+
∑
i
=
1
n
b
‍
B
i
z
-
i
]
u
-
(
t
)
,
w
(
t
)
=
[
1
+
∑
i
=
1
n
d
‍
d
i
z
-
i
]
v
(
t
)
,
(7)
y
(
t
)
=
∑
i
=
1
n
c
‍
c
i
[
γ
i
(
u
(
t
)
)
0
]
+
∑
i
=
1
n
b
‍
B
i
u
-
(
t
-
i
)
+
∑
i
=
1
n
d
‍
d
i
v
(
t
-
i
)
+
v
(
t
)
.
Note that the subscripts (Roman) s and n denote the first letters of “system” and “noise” for distinguishing the types of the unknown parameter vectors or matrices, respectively. Define the parameter matrix θ
1, the parameter vectors θ
2, θ
n, and the information vectors φ
1
(
t
), φ
2
(
t
), and φ
n
(
t
) as
(8)
θ
1
T
∶
=
[
B
1
,
B
2
,
…
,
B
n
b
]
∈
R
m
×
(
r
n
b
)
,
θ
2
∶
=
[
c
1
,
c
2
,
…
,
c
n
c
]
T
∈
R
n
c
,
θ
n
∶
=
[
d
1
,
d
2
,
…
,
d
n
d
]
T
∈
R
n
d
,
φ
1
(
t
)
∶
=
[
u
-
T
(
t
-
1
)
,
u
-
T
(
t
-
2
)
,
…
,
u
-
T
(
t
-
n
b
)
]
T
∈
R
r
n
b
,
φ
2
(
t
)
∶
=
[
η
1
(
t
)
,
η
2
(
t
)
,
…
,
η
n
c
(
t
)
]
∈
R
m
×
n
c
,
η
i
(
t
)
∶
=
[
γ
i
(
u
(
t
)
)
0
]
∈
R
m
,
φ
n
(
t
)
∶
=
[
v
(
t
-
1
)
,
v
(
t
-
2
)
,
…
,
v
(
t
-
n
d
)
]
∈
R
m
×
n
d
.
Then, we have
(9)
x
(
t
)
=
θ
1
T
φ
1
(
t
)
+
φ
2
(
t
)
θ
2
,
w
(
t
)
=
φ
n
(
t
)
θ
n
+
v
(
t
)
,
y
(
t
)
=
x
(
t
)
+
w
(
t
)
=
θ
1
T
φ
1
(
t
)
+
φ
2
(
t
)
θ
2
+
φ
n
(
t
)
θ
n
+
v
(
t
)
.
Distinguished from the hierarchical identification methods, we reparameterize the model in (7) by using the Kronecker product to get a parameter matrix ϑ and by gathering the input information vectors φ
1
(
t
) and φ
2
(
t
) and output information matrix φ
n
(
t
) into one information matrix Ψ
(
t
) to obtain an information matrix Ψ
(
t
) as follows:
(10)
ϑ
∶
=
[
θ
s
θ
n
]
∈
R
n
0
,
n
0
∶
=
m
r
n
b
+
n
c
+
n
d
,
θ
s
∶
=
[
col
[
θ
1
T
]
θ
2
]
∈
R
m
r
n
b
+
n
c
,
Ψ
(
t
)
∶
=
[
φ
s
(
t
)
,
φ
n
(
t
)
]
∈
R
m
×
n
0
,
φ
s
(
t
)
∶
=
[
φ
1
T
(
t
)
⊗
I
m
,
φ
2
(
t
)
]
∈
R
m
×
(
m
r
n
b
+
n
c
)
.
Thus, we obtain
(11)
y
(
t
)
=
Ψ
(
t
)
ϑ
+
v
(
t
)
.
Equation (11) is the identification model of the multivariable Hammerstein FIR-MA system. Defining and minimizing the cost function
(12)
J
(
ϑ
)
∶
=
∥
y
(
t
)
-
Ψ
(
t
)
ϑ
∥
2
and using the least squares search principle, we obtain the following recursive least squares algorithm [35] to obtain parameter estimates ϑ
^
(
t
):
(13)
ϑ
^
(
t
)
=
ϑ
^
(
t
-
1
)
+
L
(
t
)
[
y
(
t
)
-
Ψ
(
t
)
ϑ
^
(
t
-
1
)
]
,
(14)
L
(
t
)
=
P
(
t
-
1
)
Ψ
T
(
t
)
[
I
m
+
Ψ
(
t
)
P
(
t
-
1
)
Ψ
T
(
t
)
]
-
1
,
(15)
P
(
t
)
=
[
I
n
0
-
L
(
t
)
Ψ
(
t
)
]
P
(
t
-
1
)
,
P
(
0
)
=
p
0
I
n
0
.
Since the information matrix Ψ
(
t
) in (13) contains the unknown intermediate variables u
-
(
t
-
i
) and the unmeasurable terms v
(
t
-
i
), the recursive algorithm in (13)–(15) cannot compute the parameter estimate θ
^
(
t
). The solution here is replacing the unknown intermediate variables u
-
(
t
-
i
) and the unmeasurable terms v
(
t
-
i
) in Ψ
(
t
) with the variable estimates (or the outputs of the auxiliary model) u
-
^
(
t
-
i
) and the estimates v
^
(
t
-
i
) based on the auxiliary model identification idea. The replaced information matrices are defined as
(16)
Ψ
^
(
t
)
∶
=
[
φ
^
s
(
t
)
,
φ
^
n
(
t
)
]
∈
R
m
×
n
0
,
φ
^
s
(
t
)
∶
=
[
φ
^
1
T
(
t
)
⊗
I
m
,
φ
2
(
t
)
]
∈
R
m
×
(
m
r
n
b
+
n
c
)
,
φ
^
1
(
t
)
∶
=
[
u
-
^
T
(
t
-
1
)
,
u
-
^
T
(
t
-
2
)
,
…
,
u
-
^
T
(
t
-
n
b
)
]
T
∈
R
r
n
b
,
φ
^
n
(
t
)
∶
=
[
v
^
(
t
-
1
)
,
v
^
(
t
-
2
)
,
…
,
v
^
(
t
-
n
d
)
]
∈
R
m
×
n
d
.
Define the parameter estimation matrices
(17)
ϑ
^
(
t
)
∶
=
[
θ
^
s
(
t
)
θ
^
n
(
t
)
]
∈
R
n
0
,
θ
^
s
(
t
)
∶
=
[
col
[
θ
^
1
T
(
t
)
]
θ
^
2
(
t
)
]
∈
R
m
r
n
b
+
n
c
,
θ
^
1
T
(
t
)
∶
=
[
B
^
1
(
t
)
,
B
^
2
(
t
)
,
…
,
B
^
n
b
(
t
)
]
∈
R
m
×
(
r
n
b
)
,
θ
^
2
(
t
)
∶
=
[
c
^
1
(
t
)
,
c
^
2
(
t
)
,
…
,
c
^
n
c
(
t
)
]
T
∈
R
n
c
,
θ
^
n
(
t
)
∶
=
[
d
^
1
(
t
)
,
d
^
2
(
t
)
,
…
,
d
^
n
d
(
t
)
]
T
∈
R
n
d
.
By replacing the parameters c
i
(
i
=
1,2
,
…
,
n
c
) in (4) with c
^
i
(
t
), the output of the proposed auxiliary model u
-
^
(
t
) is given by
(18)
u
-
^
(
t
)
=
∑
i
=
1
n
c
‍
c
^
i
(
t
)
γ
i
(
u
(
t
)
)
=
c
^
1
(
t
)
γ
1
(
u
(
t
)
)
+
c
^
2
(
t
)
γ
2
(
u
(
t
)
)
+
⋯
+
c
^
n
c
(
t
)
γ
n
c
(
u
(
t
)
)
.
From (11), we obtain v
(
t
)
=
y
(
t
)
-
Ψ
(
t
)
ϑ. Replacing Ψ
(
t
) and ϑ with Ψ
^
(
t
) and ϑ
^
(
t
-
1
), the residual v
^
(
t
) can be written as v
^
(
t
)
=
y
(
t
)
-
Ψ
^
(
t
)
ϑ
^
(
t
-
1
).To summarize, we conclude the following recursive least squares algorithm for multivariable Hammerstein FIR-MA models (the MRLS algorithm for short):(19)
ϑ
^
(
t
)
=
ϑ
^
(
t
-
1
)
+
L
(
t
)
[
y
(
t
)
-
Ψ
^
(
t
)
ϑ
^
(
t
-
1
)
]
,
L
(
t
)
=
P
(
t
-
1
)
Ψ
^
T
(
t
)
[
I
m
+
Ψ
^
(
t
)
P
(
t
-
1
)
Ψ
^
T
(
t
)
]
-
1
,
P
(
t
)
=
[
I
n
0
-
L
(
t
)
Ψ
^
(
t
)
]
P
(
t
-
1
)
,
u
-
^
(
t
)
=
∑
i
=
1
n
c
‍
c
^
i
(
t
)
γ
i
(
u
(
t
)
)
,
v
^
(
t
)
=
y
(
t
)
-
Ψ
^
(
t
)
ϑ
^
(
t
-
1
)
,
Ψ
^
(
t
)
=
[
φ
^
1
T
(
t
)
⊗
I
m
,
φ
2
,
φ
^
n
(
t
)
]
,
φ
^
1
(
t
)
=
[
u
-
^
T
(
t
-
1
)
,
u
-
^
T
(
t
-
2
)
,
…
,
u
-
^
T
(
t
-
n
b
)
]
T
,
φ
2
(
t
)
=
[
η
1
(
t
)
,
η
2
(
t
)
,
…
,
η
n
c
(
t
)
]
,
η
i
(
t
)
=
[
γ
i
(
u
(
t
)
)
0
]
,
φ
^
n
(
t
)
=
[
v
^
(
t
-
1
)
,
v
^
(
t
-
2
)
,
…
,
v
^
(
t
-
n
d
)
]
,
ϑ
^
(
t
)
=
[
col
[
θ
^
1
T
(
t
)
]
θ
^
2
(
t
)
θ
^
n
(
t
)
]
,
θ
^
1
T
(
t
)
=
[
B
^
1
(
t
)
,
B
^
2
(
t
)
,
…
,
B
^
n
b
(
t
)
]
,
θ
^
2
(
t
)
=
[
c
^
1
(
t
)
,
c
^
2
(
t
)
,
…
,
c
^
n
c
(
t
)
]
T
,
θ
^
n
(
t
)
=
[
d
^
1
(
t
)
,
d
^
2
(
t
)
,
…
,
d
^
n
d
(
t
)
]
T
.


## 3. The F-MRLS Algorithm

The convergence rate of the MRLS algorithm in Section2 is slow because the noise information intermediate variables w
(
t
) contain unmeasurable time-delay noise v
(
t
-
i
). The solution here is to present a filtering based recursive least squares algorithm (the F-MRLS algorithm) for the multivariable Hammerstein models by filtering the rational function N
(
z
) and transferring the FIR-MA model in (1) into a controlled autoregressive (CAR) model. Multiplying both sides of (1) by N
-
1
(
z
) yields
(20)
N
-
1
(
z
)
y
(
t
)
=
B
(
z
)
N
-
1
(
z
)
u
-
(
t
)
+
v
(
t
)
,
or
(21)
y
f
(
t
)
=
B
(
z
)
u
-
f
(
t
)
+
v
(
t
)
,
where
(22)
u
-
f
(
t
)
∶
=
1
N
(
z
)
u
-
(
t
)
=
1
N
(
z
)
∑
i
=
1
n
c
‍
c
^
i
(
t
)
γ
i
(
u
(
t
)
)
=
c
1
U
1
(
t
)
+
c
2
U
2
(
t
)
+
⋯
+
c
n
c
U
n
c
(
t
)
,
y
f
(
t
)
∶
=
1
N
(
z
)
y
(
t
)
=
[
1
-
N
(
z
)
]
y
f
(
t
)
+
y
(
t
)
=
y
(
t
)
-
∑
i
=
1
n
d
‍
d
i
(
t
)
y
f
(
t
-
i
)
,
U
i
(
t
)
∶
=
1
N
(
z
)
γ
i
(
u
(
t
)
)
∈
R
r
,
i
=
1,2
,
…
,
n
c
,
ζ
j
(
t
)
∶
=
[
U
j
(
t
)
0
]
∈
R
m
,
j
=
1,2
,
…
,
n
c
.
Thus, (21) can be rewritten as
(23)
y
f
(
t
)
=
[
I
m
×
r
+
∑
i
=
1
n
b
‍
B
i
z
-
i
]
u
-
f
(
t
)
+
v
(
t
)
=
[
u
-
f
(
t
)
0
]
+
B
1
u
-
f
(
t
-
1
)
+
B
2
u
-
f
(
t
-
2
)
+
⋯
+
B
n
b
u
-
f
(
t
-
n
b
)
+
v
(
t
)
,
=
c
1
ζ
1
(
t
)
+
c
2
ζ
2
(
t
)
+
⋯
+
c
n
c
ζ
n
c
(
t
)
+
B
1
u
-
f
(
t
-
1
)
+
B
2
u
-
f
(
t
-
2
)
+
⋯
+
B
n
b
u
-
f
(
t
-
n
b
)
+
v
(
t
)
.
Define the filtered information matrices:
(24)
Ψ
f
(
t
)
∶
=
[
φ
f
1
T
(
t
)
⊗
I
m
,
φ
f
2
(
t
)
]
∈
R
m
×
(
m
r
n
b
+
n
c
)
,
(25)
φ
f
1
(
t
)
∶
=
[
u
-
f
T
(
t
-
1
)
,
u
-
f
T
(
t
-
2
)
,
…
,
u
-
f
T
(
t
-
n
b
)
]
T
∈
R
(
r
n
b
)
,
(26)
φ
f
2
(
t
)
∶
=
[
ζ
1
(
t
)
,
ζ
2
(
t
)
,
…
,
ζ
n
c
(
t
)
]
∈
R
m
×
n
c
.
Since the polynomial N
(
z
) is unknown and to be estimated, it is impossible to use u
-
f
(
t
) to construct φ
f
1
(
t
) in (25). Here, we adopt the principle of the MRLS algorithm in Section 2 and replace the unmeasurable variables and vectors with their estimates to derive the following algorithm.By using the parameter estimatesθ
^
1
(
t
) and θ
^
n
(
t
), the estimates of polynomials B
(
z
) and N
(
z
) at time t can be constructed as
(27)
B
^
(
t
,
z
)
∶
=
I
m
×
r
+
∑
i
=
1
n
b
‍
B
^
i
(
t
)
z
-
i
,
N
^
(
t
,
z
)
∶
=
1
+
d
^
1
(
t
)
z
-
1
+
d
^
2
(
t
)
z
-
2
+
⋯
+
d
^
n
d
(
t
)
z
-
n
d
.
w
(
t
) in model (1) can be rewritten as
(28)
w
(
t
)
=
y
(
t
)
-
x
(
t
)
=
y
(
t
)
-
[
I
m
×
r
+
∑
i
=
1
n
b
‍
B
i
(
t
)
z
-
i
]
u
-
(
t
)
=
y
(
t
)
-
φ
s
(
t
)
θ
s
.
Let w
^
(
t
) be the estimate of w
(
t
). Replacing w
(
t
), y
(
t
), φ
s
(
t
), φ
n
(
t
), θ
n, and θ
s with their estimates w
^
(
t
), y
^
(
t
), φ
^
s
(
t
), φ
^
n
(
t
), θ
^
n
(
t
), and θ
^
s
(
t
-
1
) leads to
(29)
w
^
(
t
)
=
y
(
t
)
-
φ
^
s
(
t
)
θ
^
s
(
t
-
1
)
.
Defining and minimizing the cost function
(30)
J
(
θ
n
)
∶
=
∑
j
=
1
t
‍
[
w
(
j
)
-
φ
n
(
j
)
θ
n
]
2
and using the least squares search principle, we list the recursive least squares algorithm to compute θ
^
n
(
t
):
(31)
θ
^
n
(
t
)
=
θ
^
n
(
t
-
1
)
+
L
n
(
t
)
[
w
^
(
t
)
-
φ
^
n
(
t
)
θ
^
n
(
t
-
1
)
]
,
L
n
(
t
)
=
P
n
(
t
-
1
)
φ
^
n
T
(
t
)
[
I
m
+
φ
^
n
(
t
)
P
n
(
t
-
1
)
φ
^
n
T
(
t
)
]
-
1
,
P
n
(
t
)
=
[
I
n
d
-
L
n
(
t
)
φ
^
n
(
t
)
]
P
n
(
t
-
1
)
,
P
n
(
0
)
=
p
0
I
.
Let c
^
(
t
)
∶
=
[
c
^
1
(
t
)
,
c
^
2
(
t
)
,
…
,
c
^
n
c
(
t
)
]
T
∈
R
n
c be the estimate of c at time t, filtering y
(
t
) with 1
/
N
^
(
t
,
z
) to obtain the estimate y
^
f
(
t
):
(32)
y
^
f
(
t
)
=
1
N
^
(
t
,
z
)
y
(
t
)
=
[
1
-
N
^
(
t
,
z
)
]
y
^
f
(
t
)
+
y
(
t
)
=
-
∑
i
=
1
n
d
‍
d
^
i
(
t
)
y
^
f
(
t
-
i
)
+
y
(
t
)
.
The estimate of U
j
(
t
) can be computed by
(33)
U
^
j
(
t
)
=
1
N
^
(
t
,
z
)
γ
j
(
u
(
t
)
)
=
[
1
-
N
^
(
t
,
z
)
]
U
^
j
(
t
)
+
γ
j
(
u
(
t
)
)
=
-
∑
i
=
1
n
d
‍
d
^
i
(
t
)
U
^
j
(
t
-
i
)
+
γ
j
(
u
(
t
)
)
.
Define the estimate of ζ
j
(
t
) by
(34)
ζ
^
j
(
t
)
∶
=
[
U
^
j
(
t
)
0
]
∈
R
m
,
j
=
1,2
,
…
,
n
c
,
and construct the estimate of Ψ
f
(
t
) with φ
^
f
1
(
t
) and φ
^
f
2
(
t
) as follows:
(35)
Ψ
^
f
(
t
)
=
[
φ
^
f
1
T
(
t
)
⊗
I
m
,
φ
^
f
2
(
t
)
]
,
φ
^
f
1
(
t
)
=
[
u
-
^
f
T
(
t
-
1
)
,
u
-
^
f
T
(
t
-
2
)
,
…
,
u
-
^
f
T
(
t
-
n
b
)
]
T
,
φ
^
f
2
(
t
)
=
[
ζ
^
1
(
t
)
,
ζ
^
2
(
t
)
,
…
,
ζ
^
n
c
(
t
)
]
.
The filtered model in (21) can be rewritten in a matrix form:
(36)
y
f
(
t
)
=
Ψ
f
(
t
)
θ
s
+
v
(
t
)
,
or
(37)
v
(
t
)
=
y
f
(
t
)
-
Ψ
f
(
t
)
θ
s
.
Based on the MRLS search principle, we can obtain the estimate of θ
s by the following algorithm:
(38)
θ
^
s
(
t
)
=
θ
^
s
(
t
-
1
)
+
L
f
(
t
)
[
y
^
f
(
t
)
-
Ψ
^
f
(
t
)
θ
^
s
(
t
-
1
)
]
,
L
f
(
t
)
=
P
f
(
t
-
1
)
Ψ
^
f
T
(
t
)
[
I
m
+
Ψ
^
f
(
t
)
P
f
(
t
-
1
)
Ψ
^
f
T
(
t
)
]
-
1
,
P
f
(
t
)
=
[
I
r
n
b
+
n
c
-
L
f
(
t
)
Ψ
^
f
(
t
)
]
P
f
(
t
-
1
)
,
h
h
h
h
h
h
h
h
h
h
h
i
h
h
h
h
h
P
f
(
0
)
=
p
0
I
.
The estimate u
-
^
(
t
) can be computed by
(39)
u
-
^
(
t
)
=
c
^
1
(
t
)
γ
1
(
u
(
t
)
)
+
c
^
2
(
t
)
γ
2
(
u
(
t
)
)
+
⋯
+
c
^
n
c
(
t
)
γ
n
c
(
u
(
t
)
)
.
Filter u
-
^
(
t
) by 1
/
N
^
(
t
,
z
) to obtain the estimate u
-
^
f
(
t
):
(40)
u
-
^
f
(
t
)
=
1
N
^
(
t
,
z
)
u
-
^
(
t
)
=
1
N
^
(
t
,
z
)
∑
i
=
1
n
c
‍
c
^
i
(
t
)
γ
i
(
u
(
t
)
)
=
c
^
1
(
t
)
U
^
1
(
t
)
+
c
^
2
(
t
)
U
^
2
(
t
)
+
⋯
+
c
^
n
c
(
t
)
U
^
n
c
(
t
)
.
Replacing v
(
t
), y
f
(
t
), Ψ
f
(
t
), and θ
s in (37) with their estimates v
^
(
t
), y
^
f
(
t
), Ψ
^
f
(
t
), and θ
^
s
(
t
) at time t, the noise vector can be computed by
(41)
v
^
(
t
)
=
y
^
f
(
t
)
-
Ψ
^
f
(
t
)
θ
^
s
(
t
)
.
In conclusion, we can summarize the following filtering based recursive least squares algorithm for multivariable Hammerstein models (the F-MRLS algorithm for short):
(42)
θ
^
s
(
t
)
=
θ
^
s
(
t
-
1
)
+
L
f
(
t
)
[
y
^
f
(
t
)
-
Ψ
^
f
(
t
)
θ
^
s
(
t
-
1
)
]
,
(43)
L
f
(
t
)
=
P
f
(
t
-
1
)
Ψ
^
f
T
(
t
)
I
m
+
Ψ
^
f
(
t
)
P
f
(
t
-
1
)
Ψ
^
f
T
(
t
)
,
(44)
P
f
(
t
)
=
[
I
r
n
b
+
n
c
-
L
f
(
t
)
φ
^
f
(
t
)
]
P
f
(
t
-
1
)
,
(45)
Ψ
^
f
(
t
)
=
[
φ
^
f
1
T
(
t
)
⊗
I
m
,
φ
^
f
2
]
,
(46)
φ
^
f
1
(
t
)
=
[
u
-
^
f
T
(
t
-
1
)
,
u
-
^
f
T
(
t
-
2
)
,
…
,
u
-
^
f
T
(
t
-
n
b
)
]
T
,
(47)
φ
^
f
2
(
t
)
=
[
ζ
^
1
(
t
)
,
ζ
^
2
(
t
)
,
…
,
ζ
^
n
c
(
t
)
]
,
(48)
ζ
^
j
(
t
)
=
[
U
^
j
(
t
)
0
]
,
(49)
U
^
j
(
t
)
=
-
∑
i
=
1
n
d
‍
d
^
i
(
t
)
U
^
j
(
t
-
i
)
+
γ
j
(
u
(
t
)
)
,
(50)
u
-
^
f
(
t
)
=
c
^
1
(
t
)
U
^
1
(
t
)
+
c
^
2
(
t
)
U
^
2
(
t
)
+
⋯
+
c
^
n
c
(
t
)
U
^
n
c
(
t
)
,
(51)
y
^
f
(
t
)
=
-
∑
i
=
1
n
d
‍
d
^
i
(
t
)
y
^
f
(
t
-
i
)
+
y
(
t
)
,
(52)
θ
^
n
(
t
)
=
θ
^
n
(
t
-
1
)
+
L
n
(
t
)
[
w
^
(
t
)
-
φ
^
n
(
t
)
θ
^
n
(
t
-
1
)
]
,
(53)
L
n
(
t
)
=
P
n
(
t
-
1
)
φ
^
n
T
(
t
)
[
I
m
+
φ
^
n
(
t
)
P
n
(
t
-
1
)
φ
^
n
T
(
t
)
]
-
1
,
(54)
P
n
(
t
)
=
[
I
n
d
-
L
n
(
t
)
φ
^
n
(
t
)
]
P
n
(
t
-
1
)
,
(55)
w
^
(
t
)
=
y
(
t
)
-
φ
^
s
(
t
)
θ
^
s
(
t
-
1
)
,
(56)
Ψ
^
(
t
)
=
[
φ
^
s
(
t
)
,
φ
^
n
(
t
)
]
,
(57)
φ
^
s
(
t
)
=
[
φ
^
1
T
(
t
)
⊗
I
m
,
φ
2
(
t
)
]
,
(58)
φ
^
1
(
t
)
=
[
u
-
^
T
(
t
-
1
)
,
u
-
^
T
(
t
-
2
)
,
…
,
u
-
^
T
(
t
-
n
b
)
]
T
,
(59)
φ
2
(
t
)
=
[
η
1
(
t
)
,
η
2
(
t
)
,
…
,
η
n
c
(
t
)
]
,
(60)
η
i
(
t
)
=
[
γ
i
(
u
(
t
)
)
0
]
,
(61)
φ
^
n
(
t
)
=
[
v
^
(
t
-
1
)
,
v
^
(
t
-
2
)
,
…
,
v
^
(
t
-
n
d
)
]
,
(62)
v
^
(
t
)
=
y
^
f
(
t
)
-
Ψ
^
f
(
t
)
θ
^
s
(
t
)
,
(63)
u
-
^
(
t
)
=
∑
i
=
1
n
c
‍
c
^
i
(
t
)
γ
i
(
u
(
t
)
)
,
(64)
ϑ
^
(
t
)
=
[
θ
^
s
(
t
)
θ
^
n
(
t
)
]
,
(65)
θ
^
s
(
t
)
=
[
col
[
θ
^
1
T
(
t
)
]
θ
^
2
(
t
)
]
,
(66)
θ
^
1
T
(
t
)
=
[
B
^
1
(
t
)
,
B
^
2
(
t
)
,
…
,
B
^
n
b
(
t
)
]
,
(67)
θ
^
2
(
t
)
=
[
c
^
1
(
t
)
,
c
^
2
(
t
)
,
…
,
c
^
n
c
(
t
)
]
T
,
(68)
θ
^
n
(
t
)
=
[
d
^
1
(
t
)
,
d
^
2
(
t
)
,
…
,
d
^
n
d
(
t
)
]
T
.The steps involved in the F-MRLS algorithm for multivariable Hammerstein systems are listed in the following.(1)
To initialize, lett
=
1, set the initial values of the parameter estimation variables and covariance matrices as follows: θ
^
s
(
i
)
=
1
m
r
n
b
+
n
c
/
p
0, θ
^
n
(
i
)
=
1
n
d
/
p
0, y
^
f
(
i
)
=
1
m
/
p
0, u
-
^
f
(
i
)
=
1
r
/
p
0, u
-
^
(
i
)
=
1
r
/
p
0, v
^
(
i
)
=
1
m
/
p
0 for i
⩽
0, U
^
j
(
i
)
=
1
r
/
p
0 for i
⩽
0 and j
=
1,2
,
…
,
n
c, P
f
(
0
)
=
p
0
I
r
n
b
+
n
c, P
n
(
0
)
=
p
0
I
n
d, p
0
=
1
0
6, and give the basis functions γ
j
(
·
).
(2)
Collect the input-output datau
(
t
) and y
(
t
), and construct η
i
(
t
) by (60), and φ
^
1
(
t
), φ
2
(
t
) by (58), (59). Form the information vectors φ
^
s
(
t
) by (57), φ
^
n
(
t
) by (61), and Ψ
^
(
t
) by (56), respectively.
(3)
Computew
^
(
t
) by (55), the gain vector L
n
(
t
) by (53), and the covariance matrix P
n
(
t
) by (54), respectively. Update the parameter estimate θ
^
n
(
t
) by (52).
(4)
Computey
^
f
(
t
) by (51), U
^
j
(
t
) by (49). Construct ζ
^
j
(
t
), φ
^
f
2
(
t
), and φ
^
f
1
(
t
) by (48), (47), and (46), respectively. Compute Ψ
^
f
(
t
) by (45).
(5)
Compute the gain vectorL
f
(
t
) by (43) and the covariance matrix P
f
(
t
) by (44). Update the parameter estimate θ
^
s
(
t
) by (42).
(6)
Construct the parameter vectorsθ
^
1
(
t
), θ
^
2
(
t
), and θ
^
n
(
t
) by (66), (67), and (68). Form θ
^
s
(
t
) and ϑ
^
(
t
) by (65) and (64). Compute u
-
^
f
(
t
) by (50), v
^
(
t
) by (62), and u
-
^
(
t
) by (63).
(7)
Increaset by 1 and go to step 2.


## 4. Examples

Example 1.
Consider the following 2-input 2-output Hammerstein FIR-MA system:(69)
[
y
1
(
t
)
y
2
(
t
)
]
=
[
u
-
1
(
t
)
u
-
2
(
t
)
]
+
B
[
u
-
1
(
t
-
1
)
u
-
2
(
t
-
1
)
]
+
[
v
1
(
t
)
v
2
(
t
)
]
+
d
1
[
v
1
(
t
-
1
)
v
2
(
t
-
1
)
]
,
where
(70)
B
(
z
)
=
[
0.13
0.25
-
1.21
0.13
]
z
-
1
,
D
(
z
)
=
1
+
d
1
z
-
1
=
1
+
0.68
z
-
1
,
u
-
(
t
)
=
c
1
u
(
t
)
+
c
2
u
2
(
t
)
,
c
=
[
c
1
,
c
2
]
T
=
[
1.00
,
-
0.14
]
T
,
ϑ
=
[
0.13
,
-
1.21,0.25,0.13
,
-
0.14,0.68
]
T
.
In simulation, the inputs {
u
1
(
t
)
} and {
u
2
(
t
)
} are taken as two uncorrelated persistent excitation signal sequences with zero mean and unit variances and {
v
1
(
t
)
} and {
v
2
(
t
)
} as two white noise sequences with zero mean and variances σ
1
2
=
σ
2
2
=
σ
2
(
σ
2
=
0.5
0
2
,
σ
2
=
1.0
0
2
). Applying the MRLS algorithm in (19) and the F-MRLS algorithm in (42)–(68) to estimate the parameters of this multivariable Hammerstein system, the F-MRLS parameter estimates and their estimation errors are shown in Table 1, the comparison between the F-MRLS algorithm and the MRLS algorithm in the estimation error δ
∶
=
∥
ϑ
^
(
t
)
-
ϑ
∥
/
∥
ϑ
∥ versus t is shown in Table 2, and the estimation errors δ
∶
=
∥
ϑ
^
(
t
)
-
ϑ
∥
/
∥
ϑ
∥ versus t are shown in Figures 1 and 2.Table 1
The F-MRLS estimates and errors in Example1 (σ
2
=
0.5
0
2 and σ
2
=
1.0
0
2).
σ
2
t
B
(
1,1
)
B
(
2,1
)
B
(
1,2
)
B
(
2,2
)
c
2
d
1
δ (%)
0.5
0
2
100
0.10919
−1.60279
0.39771
0.12515
−0.16795
0.46319
33.14292
200
0.08155
−1.44298
0.31145
0.09590
−0.13875
0.51078
21.01533
500
0.10144
−1.33934
0.25923
0.08705
−0.15442
0.60335
11.18663
1000
0.12452
−1.30385
0.26703
0.11599
−0.14796
0.65282
7.04142
2000
0.13930
−1.27241
0.26711
0.12703
−0.14550
0.69136
4.66378
3000
0.14008
−1.25781
0.26567
0.11928
−0.14361
0.70371
4.03373
1.0
0
2
100
0.09530
−1.72159
0.52726
0.13387
−0.15869
0.52551
42.21841
200
0.03631
−1.52223
0.36553
0.07113
−0.11792
0.58022
25.56872
500
0.07400
−1.40400
0.25492
0.04715
−0.16037
0.64763
15.50868
1000
0.12015
−1.36311
0.27803
0.10606
−0.15335
0.68001
11.08080
2000
0.14849
−1.31457
0.28056
0.12704
−0.15048
0.70624
7.98334
3000
0.15005
−1.29138
0.27909
0.10998
−0.14680
0.71335
6.79480
True values
0.13000
−1.21000
0.25000
0.13000
−0.14000
0.68000Table 2
The MRLS and F-MRLS estimates and errors in Example1 (σ
2
=
0.5
0
2).
Algorithms
t
B
(
1,1
)
B
(
2,1
)
B
(
1,2
)
B
(
2,2
)
c
2
d
1
δ (%)
MRLS
100
−0.01946
−1.84785
0.44644
0.04712
−0.16528
0.39802
52.11950
200
−0.03340
−1.60645
0.31299
0.05091
−0.14456
0.46017
34.45327
500
0.03788
−1.39775
0.28600
0.03749
−0.14940
0.56178
18.20046
1000
0.08074
−1.34458
0.28178
0.07949
−0.14078
0.61722
11.71737
2000
0.11167
−1.28751
0.27285
0.10106
−0.13785
0.66411
6.24292
3000
0.12600
−1.27455
0.26866
0.10656
−0.13729
0.68367
4.99772
F-MRLS
100
0.10919
−1.60279
0.39771
0.12515
−0.16795
0.46319
33.14292
200
0.08155
−1.44298
0.31145
0.09590
−0.13875
0.51078
21.01533
500
0.10144
−1.33934
0.25923
0.08705
−0.15442
0.60335
11.18663
1000
0.12452
−1.30385
0.26703
0.11599
−0.14796
0.65282
7.04142
2000
0.13930
−1.27241
0.26711
0.12703
−0.14550
0.69136
4.66378
3000
0.14008
−1.25781
0.26567
0.11928
−0.14361
0.70371
4.03373
True values
0.13000
−1.21000
0.25000
0.13000
−0.14000
0.68000Figure 1
The F-RLS estimation errorsδ versus t in Example 1 (σ
2
=
0.5
0
2 and σ
2
=
1.0
0
2).Figure 2
The parameter estimation errorsδ versus t in Example 1 (σ
2
=
0.5
0
2).Example 2.
Consider the following 2-input 2-output Hammerstein FIR-MA system:(71)
[
y
1
(
t
)
y
2
(
t
)
]
=
[
u
-
1
(
t
)
u
-
2
(
t
)
]
+
B
[
u
-
1
(
t
-
1
)
u
-
2
(
t
-
1
)
]
+
[
v
1
(
t
)
v
2
(
t
)
]
+
d
1
[
v
1
(
t
-
1
)
v
2
(
t
-
1
)
]
,
where
(72)
B
(
z
)
=
[
0.14
0.25
-
0.15
0.125
]
z
-
1
,
D
(
z
)
=
1
+
d
1
z
-
1
=
1
+
0.50
z
-
1
,
u
-
(
t
)
=
c
1
u
(
t
)
+
c
2
u
2
(
t
)
+
c
3
u
3
(
t
)
,
c
=
[
c
1
,
c
2
,
c
3
]
T
=
[
1.00
,
-
0.19,1.19
]
T
,
ϑ
=
[
0.14
,
-
0.15,0.25,0.125
,
-
0.19,1.19,0.50
]
T
.
In simulation, the inputs {
u
1
(
t
)
} and {
u
2
(
t
)
} and noise data {
v
1
(
t
)
} and {
v
2
(
t
)
} are settled in the same way as in Example 1, where variances σ
1
2
=
σ
2
2
=
σ
2
(
σ
2
=
0.5
0
2). Applying the MRLS algorithm and the F-MRLS algorithm in (42)–(68), the F-MRLS parameter estimates and their estimation errors are shown in Table 3, the comparison between the F-MRLS algorithm and the MRLS algorithm in the estimation error δ
∶
=
∥
ϑ
^
(
t
)
-
ϑ
∥
/
∥
ϑ
∥ versus t is shown in Table 4, and the estimation errors δ
∶
=
∥
ϑ
^
(
t
)
-
ϑ
∥
/
∥
ϑ
∥ versus t are shown in Figures 3 and 4.Table 3
The F-RLS estimates and errors in Example2 (σ
2
=
0.5
0
2 and σ
2
=
1.0
0
2).
σ
2
t
B
(
1,1
)
B
(
2,1
)
B
(
1,2
)
B
(
2,2
)
c
2
c
3
d
1
δ (%)
0.5
0
2
100
0.12373
−1.02307
0.27283
0.10358
−0.39300
1.25093
0.32512
16.33820
200
0.12306
−1.03981
0.24977
0.11748
−0.29306
1.19925
0.35407
10.57344
500
0.13109
−1.05853
0.24655
0.11563
−0.25013
1.19098
0.40948
6.44687
1000
0.13728
−1.06400
0.25207
0.12027
−0.22484
1.21159
0.45612
3.63495
2000
0.14372
−1.06203
0.25315
0.12346
−0.20822
1.20110
0.49970
1.46876
3000
0.14389
−1.05999
0.25321
0.12156
−0.20281
1.19897
0.51458
1.43078
1.0
0
2
100
0.10227
−0.88895
0.29620
0.02428
−0.45977
1.39972
0.43173
23.55652
200
0.10022
−0.94211
0.26126
0.06755
−0.30135
1.31216
0.41341
13.31668
500
0.11893
−1.01753
0.24656
0.08815
−0.27093
1.24343
0.43981
7.40054
1000
0.13195
−1.05047
0.25483
0.10514
−0.24217
1.26386
0.47397
5.67199
2000
0.14551
−1.05983
0.25634
0.11672
−0.21751
1.22939
0.51038
3.02299
3000
0.14632
−1.06041
0.25650
0.11459
−0.20918
1.22061
0.52227
2.68962
True values
0.14000
−1.05000
0.25000
0.12500
−0.19000
1.19000
0.50000Table 4
The comparison of parameter estimates and errors in Example2 (σ
2
=
0.5
0
2).
Algorithms
t
B
(
1,1
)
B
(
2,1
)
B
(
1,2
)
B
(
2,2
)
c
2
c
3
d
1
δ (%)
MRLS
100
0.10379
−1.49022
0.30877
0.11454
−0.37596
0.96745
0.25108
34.46417
200
0.08551
−1.31358
0.29070
0.10159
−0.30829
1.05130
0.24601
24.37082
500
0.10950
−1.17837
0.27162
0.09651
−0.24919
1.13023
0.30548
14.80320
1000
0.12218
−1.13050
0.26705
0.10410
−0.21956
1.17878
0.36159
9.76634
2000
0.13294
−1.09517
0.26234
0.11047
−0.20567
1.18593
0.42481
5.36992
3000
0.13715
−1.08439
0.25990
0.11297
−0.20141
1.18517
0.45571
3.49648
F-MRLS
100
0.12373
−1.02307
0.27283
0.10358
−0.39300
1.25093
0.32512
16.33820
200
0.12306
−1.03981
0.24977
0.11748
−0.29306
1.19925
0.35407
10.57344
500
0.13109
−1.05853
0.24655
0.11563
−0.25013
1.19098
0.40948
6.44687
1000
0.13728
−1.06400
0.25207
0.12027
−0.22484
1.21159
0.45612
3.63495
2000
0.14372
−1.06203
0.25315
0.12346
−0.20822
1.20110
0.49970
1.46876
3000
0.14389
−1.05999
0.25321
0.12156
−0.20281
1.19897
0.51458
1.43078
True values
0.14000
−1.05000
0.25000
0.12500
−0.19000
1.19000
0.50000Figure 3
The F-RLS estimation errorsδ versus t in Example 2 (σ
2
=
0.5
0
2 and σ
2
=
1.0
0
2).Figure 4
The parameter estimation errorsδ versus t in Example 2 (σ
2
=
0.5
0
2).To illustrate the advantages of the proposed algorithm, the numbers of multiplications and additions for each step of the F-MRLS algorithm and the MRLS algorithm are listed in Table5.Table 5
Comparison of the computational efficiency of the F-MRLS and MRLS algorithms.
Algorithms
Number of multiplications
Number of additions
MRLS
2
n
2
+
5
n
2
n
2
+
2
n
F-MRLS
2
n
b
2
+
2
n
c
2
+
2
n
d
2
+
5
n
2
n
b
2
+
2
n
c
2
+
2
n
d
2
+
2
n
+
1From Tables1–5, Figures 1–4, we can draw the following conclusions.(i)
The parameter estimation errors are getting smaller witht increasing, which proves that the proposed algorithms are effective.
(ii)
The F-MRLS algorithm is more accurate than the MRLS algorithm, which means the proposed F-MRLS algorithm has a better performance compared with the MRLS algorithm.
(iii)
The parameter estimates given by the F-MRLS algorithm have faster convergence than those given by the MRLS algorithm.


## 5. Conclusions

This paper presents a data filtering based recursive least squares algorithm for MIMO nonlinear FIR-MA systems. The simulation results show that the proposed data filtering based recursive least squares algorithm is more accurate and reduces computational burden compared with the recursive least squares algorithm.



---

*Source: 232848-2014-10-16.xml*