# Recursive Subspace Identification of AUV Dynamic Model under General Noise Assumption

Authors: Zheping Yan; Di Wu; Jiajia Zhou; Lichao Hao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A recursive subspace identification algorithm for autonomous underwater vehicles (AUVs) is proposed in this paper. Due to the advantages at handling nonlinearities and couplings, the AUV model investigated here is for the first time constructed as a Hammerstein model with nonlinear feedback in the linear part. To better take the environment and sensor noises into consideration, the identification problem is concerned as an errors-in-variables (EIV) one which means that the identification procedure is under general noise assumption. In order to make the algorithm recursively, propagator method (PM) based subspace approach is extended into EIV framework to form the recursive identification method called PM-EIV algorithm. With several identification experiments carried out by the AUV simulation platform, the proposed algorithm demonstrates its effectiveness and feasibility.


---

## Body

## 1. Introduction

In recent years, autonomous underwater vehicles (AUVs) have attracted increasing attentions due to their remarkable features such as high agility, excellent convenience and low cost in applications of underwater explorations and developments. However, contradictions lay between more and more complicated missions for AUV and the control and navigation systems that are not accurate enough. System identification methods have provided an alternative way other than traditional expensive instruments dependent approaches to improve the abilities of autonomous systems in various aspects [1–5]. As a result, a variety of researches have been put forward to identify the ordinary differential model of AUVs for model based control and navigation. Rentschler and coworkers [1] have demonstrated an iterative procedure to revise the model and controller of Odyssey III AUV to obtain better flight performances. Nonlinear observers based identification algorithm with sliding mode observer and EKF is also proposed for designation of nonlinear controller in [2]. For a more robust navigation system in case of sensor fault, Hegrenaes and Hallingstad [3] have used least squares algorithm to estimate both sea current disturbances and model parameters of the HUGIN 4500 AUV. However, due to the complexity of AUV system, many nonlinearities and coupled terms exist in ordinary differential equations that make the identification of whole pack of hydrodynamic coefficients quite complicated and time consuming. As a consequence, AUV differential model usually need to be simplified through eliminations of nonlinear and coupled terms before identification process. For example, in the research of Tiano et al. [6], a set of decoupled AUV models concerning different degree of freedom were set up and the yaw dynamics of the Hammerhead AUV was identified according to observer Kalman filter identification (OKID) method. However, nonlinearities and couplings are two significant features being researched in the area of AUV, so as it is said in [6], construction of MIMO coupled AUV model is a rather important and challenging modeling issue.In this paper, compared with the traditional differential equation used in AUV modeling, a Hammerstein model which consists of a static nonlinear part and a dynamic linear one is adopted in order to deal with nonlinear and linear property of AUV system separately. Due to the particular characteristics of AUV system, the Hammerstein model has to be modified with a static nonlinear feedback part added on the linear part. As illustrated in Figure1, one remarkable benefit of Hammerstein system is that it possesses advantages of linear MIMO system and can approximate nonlinear and coupled terms of AUV to a large extent at the same time. To the best of the authors’ acknowledge, this is the first time that Hammerstein model is applied in the area of AUV modeling. Because of the conveniences brought by Hammerstein system, extensive attentions has been paid to obtain system parameters from input/output data. Cai and Bai [7] proposed a method to make the identification of parametric Hammerstein system a linear problem through regarding the average squared error cost function as the inner product between the true but unknown parameter vector and its estimations. But this method was discussed under the assumption that Hammerstein system only consists of single input and single output. As to MIMO Hammerstein system identification, a nonparametric algorithm based on stochastic approximation approach is proposed in [8]. However, nonlinear MIMO Hammerstein identification problem concerned in this paper is still challenging.Figure 1
Hammerstein AUV Model.One group of widely studied MIMO system identification algorithms for Hammerstein system is subspace identification methods [9], which mainly include three different branches: numerical algorithms for state-space subspace system identification (N4SID), MIMO output-error state-space model identification (MOESP) and canonical variate analysis (CVA). Compared with other identification algorithms, subspace identification methods are more attractive due to several advantages [10]. For example, subspace identification methods can circumvent the complicated parameterization procedure for MIMO system of prediction error methods [11]. What is more, subspace identification methods do not require nonlinear searches in the parameter space based on computational tools such as QR factorization and singular value decomposition (SVD) [12]. Since the Hammerstein AUV model is a parametric one with multiinputs and multioutputs, MIMO MOESP algorithm is adopted as the theoretical basis for further investigation in this paper.In addition, instead of regarding the identification procedure in ideal situations, practical engineering circumstances need to be taken into consideration. A widely studied one is that the general noise assumption has to be made which means “the measured input is corrupted by a white measurement noise while the measured output is corrupted by the sum of a white measurement noise and a term due to a white process noise” [13]. Another practical situation which need to be considered is that different oceanic environment may lead to different hydrodynamic coefficients in which case off-line identification results of Hammerstein AUV model will bring errors to control and navigation systems [14]. Besides, it is often the case that the structure of AUV usually has to be modified mildly in order to fulfill various tasks in practice. Therefore, recursive identification methods which can adjust model parameters online become rather significant and attractive in applications of such as adaptive control, model-aided navigation and so on.As a result, in order to fulfill the situations concerned above, the identification of Hammerstein AUV model is regarded as an errors-in-variables (EIV) problem in this paper and we are going to make the MOESP method recursively under EIV framework so, as to carry out recursive identification algorithm for the Hammerstein AUV model. One major obstacle in recursive subspace identification is the increasing computational complexity of SVD [15, 16]. In previous studies, projection approximation subspace tracking (PAST) algorithm designed by Yang [17] was widely used to update the SVD. However, approximation in the algorithm will bring slight difference between identified model and the original one. In addition, IV-PAST algorithm [18] and gradient type subspace tracking method [19] are also developed to estimate the signal subspace. In [20], an instrumental variable propagator method (IVPM) based recursive identification algorithm was proposed. Nevertheless, the paper only studied the method within past input (PI)/past output (PO) MOESP framework. So in this paper, combining with Hammerstein AUV identification problem described above, IVPM method is extended into EIV framework and the PM-EIV algorithm for recursive subspace identification of errors-in-variables problem is derived. Compared with previous algorithms mentioned, the PM-EIV algorithm is more suitable to handle the MIMO AUV Hammerstein model identification problems under errors-in-variables framework recursively. As a matter of fact, since the Hammerstein model constructed can be viewed as a generalized one for mechanical systems, the proposed PM-EIV can be extended to identification of other systems as well.The remainder of this paper is organized as follows. In Section2, a Hammerstein AUV model with nonlinear feedback is formed with proper transformation from an ordinary differential one. Then the linearization process of nonlinear part is presented. In Section 3, the MOESP identification method is described with no consideration about the system noise. After that, a recursive subspace method called PM-EIV is derived under the general noise assumption. Identification and validation of the Hammerstein AUV model are presented in Section 4. At last, conclusions are made in Section 5.Some notations used in this paper are followed. The superscript(
·
)
T denotes the transposition operator. E
(
·
) is the expectation operator. R
m
×
n represents the set of m
×
n real matrices. M
⊥ is the orthogonal complement of M.


## 2. AUV Modeling

In this section, a Hammerstein AUV model is formed after introducing the ordinary differential one. Then in order to make the Hammerstein AUV model suitable for subspace identification method, linearization of the nonlinear part is brought out.


### 2.1. Hammerstein AUV Model

According to [21], a coupled nonlinear ordinary differential model for AUV based on Newtonian mechanics can be described as the following equation:
(1)
M
ν
˙
+
C
(
ν
)
ν
+
D
(
ν
)
ν
=
τ
,
where ν
∈
R
6 represents the state vector of AUV, M
∈
R
6
×
6 consists of inertia matrix and add mass matrix, D
(
ν
)
∈
R
6
×
6 is the linear and quadratic damping matrix, C
(
ν
)
∈
R
6
×
6 is the coriolis and centripetal matrix, τ
∈
R
6 represents the forces and moments acted on the vehicle.The coupled and nonlinear terms in matrixD
(
ν
) and C
(
ν
) make the identification process quite complicated and time consuming. To obtain Hammerstein model of AUV, simplifications and transformations have to be made. First, some coupled terms with little influences are eliminated. Second, the nonlinear terms are separated from the system and the system is divided into static nonlinear input part and nonlinear feedback part. Then, remaining nonlinear and coupled terms constitute the nonlinear feedback part. At last, the remaining part of the model can be described as a linear MIMO state-space model. After those steps, a Hammerstein AUV model with nonlinear feedback part can be formed as in Figure 1. d
(
k
), υ
(
k
), and ω
(
k
) are process noise, input measurement noise, and output measurement noise respectively.According to the system structure above, define the nonlinear input function asf
(
·
), nonlinear feedback as g
(
·
), so discrete time MIMO state-space equations can be constructed as below. 								(2)
x
(
k
+
1
)
=
A
x
(
k
)
+
B
1
f
(
u
(
k
)
)
+
B
2
g
(
y
(
k
)
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
D
1
f
(
u
(
k
)
)
+
D
2
g
(
y
(
k
)
)
,
y
~
(
k
)
=
y
(
k
)
+
ω
(
k
)
,
u
~
(
k
)
=
u
(
k
)
+
υ
(
k
)
,
where x
(
k
)
∈
R
6 represent the system states, u
(
k
)
∈
R
3 are the system inputs including thruster revolution, rudder deflection, and elevator deflection, y
(
k
)
∈
R
6 is output vector, consists of u, v, w, p, q, r that represents the surge velocity, sway velocity, heave velocity, roll rate, pitch rate, yaw rate respectively. A, B
1, B
2, C, D
1, and D
2 are the corresponding linear subsystem matrices with appropriate dimensions. u
~
(
k
), y
~
(
k
) are corrupted system inputs and outputs.Many identification methods have focused on handling output measurement noise and process noise. In fact, input measurement noise is inevitable in any engineering processes including practical identification experiments. So in this paper,d
(
k
), υ
(
k
), and ω
(
k
) are considered under the general noise assumption.


### 2.2. Linearization of AUV Model

According to the Section2.1, AUV model can be described as Hammerstein one with nonlinear feedback. In order to identify the parameters of corresponding matrices, nonlinearity in the equations needs to be linearly parameterized so that recursive subspace identification methods can be applied. One traditional approach to approximate nonlinearity of the system is linear combination of basic functions. In Lovera [12], Tchebiceff polynomials are chosen to linearize nonlinearities. However, as the system studied in this paper contains nonlinear feedback part, Tchebiceff polynomials approximation of g
(
·
) will influence the identifiability of system matrix. So in this paper, truncated Fourier series described in Luo and Leonessa [22], also known as trigonometric polynomials are adopted to linearize the nonlinear functions f
(
·
) and g
(
·
). Define the basic function as follows:
(3)
φ
k
(
x
)
=
1
,
k
=
0
,
φ
k
(
x
)
=
[
cos
⁡
(
k
π
(
x
-
x
m
)
x
d
)
sin
(
k
π
(
x
-
x
m
)
x
d
)
]
T
,
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
k
≥
1
.Then nonlinear function can be approximated by following equation:(4)
F
(
x
)
=
w
o
+
∑
i
=
1
N
w
i
φ
i
(
x
)
,
where x
m
=
(
x
max
⁡
+
x
min
⁡
)
/
2,  x
d
=
(
x
max
⁡
-
x
min
⁡
)
/
2,  and w
i
=
[
w
i
cos
⁡
w
i
sin
]. Since the Hammerstein AUV model is MIMO, define Φ
(
x
)
=
[
φ
0
T
(
x
)
φ
1
T
(
x
)
⋯
φ
N
T
(
x
)
]
T,  W
=
[
w
0
w
1
⋯
w
N
]; then f
(
·
), g
(
·
) can be expressed as:
(5)
f
(
u
)
=
[
W
u
1
T
W
u
2
T
⋯
W
u
r
T
]
T
·
[
Φ
T
(
u
1
)
Φ
T
(
u
2
)
⋯
Φ
T
(
u
m
)
]
T
,
g
(
y
)
=
[
W
y
1
T
W
y
2
T
⋯
W
y
s
T
]
T
·
[
Φ
T
(
y
1
)
Φ
T
(
y
2
)
⋯
Φ
T
(
y
l
)
]
T
,
where m
=
3, l
=
6, and corresponding coefficient vector W
u
i
∈
R
1
×
3
(
2
N
+
1
) and W
y
i
∈
R
1
×
6
(
2
N
+
1
). A simplified expression for nonlinear part can be presented as follows:
(6)
f
(
u
)
=
K
ξ
(
u
)
,
g
(
y
)
=
P
ζ
(
y
)
,
with definitions that ξ
(
u
)
=
Δ
[
Φ
T
(
u
1
)
Φ
T
(
u
2
)
⋯
Φ
T
(
u
m
)
]
T, K
=
Δ
[
W
u
1
T
W
u
2
T
⋯
W
u
r
T
]
T. The equation for g
(
y
) can also be formed in a similar principle. Then a new state-space model with no nonlinearity can be described as follow:
(7)
x
(
k
+
1
)
=
A
x
(
k
)
+
M
1
ξ
(
u
(
k
)
)
+
M
2
ζ
(
y
(
k
)
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
N
1
ξ
(
u
(
k
)
)
+
N
2
ζ
(
y
(
k
)
)
,
where M
1
=
B
1
K
∈
R
n
×
3
(
2
N
+
1
),  M
2
=
B
2
P
∈
R
n
×
6
(
2
N
+
1
),  N
1
=
D
1
K
∈
R
l
×
3
(
2
N
+
1
), and  N
2
=
D
2
P
∈
R
l
×
6
(
2
N
+
1
) represent new coefficients matrices, ξ
(
u
(
k
)
)
∈
R
3
(
2
N
+
1
)
×
1,  ζ
(
y
(
k
)
)
∈
R
6
(
2
N
+
1
)
×
1, respectively represent the input vector and the feedback vector. Further, define M
=
[
M
1
M
2
],  N
=
[
N
1
N
2
],  e
(
k
)
=
[
ξ
T
(
u
(
k
)
)
ζ
T
(
y
(
k
)
)
]
T as the system matrix and input vector and the Hammerstein AUV model can be presented in classical form:
(8)
x
(
k
+
1
)
=
A
x
(
k
)
+
M
e
(
k
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
N
e
(
k
)
.By now, a Hammerstein AUV model is established and has been linearized based on trigonometric polynomials with a description as in (8) being obtained. In the following section, a recursive subspace method to identify the system under EIV frame will be derived.


## 2.1. Hammerstein AUV Model

According to [21], a coupled nonlinear ordinary differential model for AUV based on Newtonian mechanics can be described as the following equation:
(1)
M
ν
˙
+
C
(
ν
)
ν
+
D
(
ν
)
ν
=
τ
,
where ν
∈
R
6 represents the state vector of AUV, M
∈
R
6
×
6 consists of inertia matrix and add mass matrix, D
(
ν
)
∈
R
6
×
6 is the linear and quadratic damping matrix, C
(
ν
)
∈
R
6
×
6 is the coriolis and centripetal matrix, τ
∈
R
6 represents the forces and moments acted on the vehicle.The coupled and nonlinear terms in matrixD
(
ν
) and C
(
ν
) make the identification process quite complicated and time consuming. To obtain Hammerstein model of AUV, simplifications and transformations have to be made. First, some coupled terms with little influences are eliminated. Second, the nonlinear terms are separated from the system and the system is divided into static nonlinear input part and nonlinear feedback part. Then, remaining nonlinear and coupled terms constitute the nonlinear feedback part. At last, the remaining part of the model can be described as a linear MIMO state-space model. After those steps, a Hammerstein AUV model with nonlinear feedback part can be formed as in Figure 1. d
(
k
), υ
(
k
), and ω
(
k
) are process noise, input measurement noise, and output measurement noise respectively.According to the system structure above, define the nonlinear input function asf
(
·
), nonlinear feedback as g
(
·
), so discrete time MIMO state-space equations can be constructed as below. 								(2)
x
(
k
+
1
)
=
A
x
(
k
)
+
B
1
f
(
u
(
k
)
)
+
B
2
g
(
y
(
k
)
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
D
1
f
(
u
(
k
)
)
+
D
2
g
(
y
(
k
)
)
,
y
~
(
k
)
=
y
(
k
)
+
ω
(
k
)
,
u
~
(
k
)
=
u
(
k
)
+
υ
(
k
)
,
where x
(
k
)
∈
R
6 represent the system states, u
(
k
)
∈
R
3 are the system inputs including thruster revolution, rudder deflection, and elevator deflection, y
(
k
)
∈
R
6 is output vector, consists of u, v, w, p, q, r that represents the surge velocity, sway velocity, heave velocity, roll rate, pitch rate, yaw rate respectively. A, B
1, B
2, C, D
1, and D
2 are the corresponding linear subsystem matrices with appropriate dimensions. u
~
(
k
), y
~
(
k
) are corrupted system inputs and outputs.Many identification methods have focused on handling output measurement noise and process noise. In fact, input measurement noise is inevitable in any engineering processes including practical identification experiments. So in this paper,d
(
k
), υ
(
k
), and ω
(
k
) are considered under the general noise assumption.


## 2.2. Linearization of AUV Model

According to the Section2.1, AUV model can be described as Hammerstein one with nonlinear feedback. In order to identify the parameters of corresponding matrices, nonlinearity in the equations needs to be linearly parameterized so that recursive subspace identification methods can be applied. One traditional approach to approximate nonlinearity of the system is linear combination of basic functions. In Lovera [12], Tchebiceff polynomials are chosen to linearize nonlinearities. However, as the system studied in this paper contains nonlinear feedback part, Tchebiceff polynomials approximation of g
(
·
) will influence the identifiability of system matrix. So in this paper, truncated Fourier series described in Luo and Leonessa [22], also known as trigonometric polynomials are adopted to linearize the nonlinear functions f
(
·
) and g
(
·
). Define the basic function as follows:
(3)
φ
k
(
x
)
=
1
,
k
=
0
,
φ
k
(
x
)
=
[
cos
⁡
(
k
π
(
x
-
x
m
)
x
d
)
sin
(
k
π
(
x
-
x
m
)
x
d
)
]
T
,
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
o
k
≥
1
.Then nonlinear function can be approximated by following equation:(4)
F
(
x
)
=
w
o
+
∑
i
=
1
N
w
i
φ
i
(
x
)
,
where x
m
=
(
x
max
⁡
+
x
min
⁡
)
/
2,  x
d
=
(
x
max
⁡
-
x
min
⁡
)
/
2,  and w
i
=
[
w
i
cos
⁡
w
i
sin
]. Since the Hammerstein AUV model is MIMO, define Φ
(
x
)
=
[
φ
0
T
(
x
)
φ
1
T
(
x
)
⋯
φ
N
T
(
x
)
]
T,  W
=
[
w
0
w
1
⋯
w
N
]; then f
(
·
), g
(
·
) can be expressed as:
(5)
f
(
u
)
=
[
W
u
1
T
W
u
2
T
⋯
W
u
r
T
]
T
·
[
Φ
T
(
u
1
)
Φ
T
(
u
2
)
⋯
Φ
T
(
u
m
)
]
T
,
g
(
y
)
=
[
W
y
1
T
W
y
2
T
⋯
W
y
s
T
]
T
·
[
Φ
T
(
y
1
)
Φ
T
(
y
2
)
⋯
Φ
T
(
y
l
)
]
T
,
where m
=
3, l
=
6, and corresponding coefficient vector W
u
i
∈
R
1
×
3
(
2
N
+
1
) and W
y
i
∈
R
1
×
6
(
2
N
+
1
). A simplified expression for nonlinear part can be presented as follows:
(6)
f
(
u
)
=
K
ξ
(
u
)
,
g
(
y
)
=
P
ζ
(
y
)
,
with definitions that ξ
(
u
)
=
Δ
[
Φ
T
(
u
1
)
Φ
T
(
u
2
)
⋯
Φ
T
(
u
m
)
]
T, K
=
Δ
[
W
u
1
T
W
u
2
T
⋯
W
u
r
T
]
T. The equation for g
(
y
) can also be formed in a similar principle. Then a new state-space model with no nonlinearity can be described as follow:
(7)
x
(
k
+
1
)
=
A
x
(
k
)
+
M
1
ξ
(
u
(
k
)
)
+
M
2
ζ
(
y
(
k
)
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
N
1
ξ
(
u
(
k
)
)
+
N
2
ζ
(
y
(
k
)
)
,
where M
1
=
B
1
K
∈
R
n
×
3
(
2
N
+
1
),  M
2
=
B
2
P
∈
R
n
×
6
(
2
N
+
1
),  N
1
=
D
1
K
∈
R
l
×
3
(
2
N
+
1
), and  N
2
=
D
2
P
∈
R
l
×
6
(
2
N
+
1
) represent new coefficients matrices, ξ
(
u
(
k
)
)
∈
R
3
(
2
N
+
1
)
×
1,  ζ
(
y
(
k
)
)
∈
R
6
(
2
N
+
1
)
×
1, respectively represent the input vector and the feedback vector. Further, define M
=
[
M
1
M
2
],  N
=
[
N
1
N
2
],  e
(
k
)
=
[
ξ
T
(
u
(
k
)
)
ζ
T
(
y
(
k
)
)
]
T as the system matrix and input vector and the Hammerstein AUV model can be presented in classical form:
(8)
x
(
k
+
1
)
=
A
x
(
k
)
+
M
e
(
k
)
+
d
(
k
)
,
y
(
k
)
=
C
x
(
k
)
+
N
e
(
k
)
.By now, a Hammerstein AUV model is established and has been linearized based on trigonometric polynomials with a description as in (8) being obtained. In the following section, a recursive subspace method to identify the system under EIV frame will be derived.


## 3. Recursive Subspace Identification Method

In this section, the PM-EIV subspace identification method for the Hammerstein AUV model will be proposed. To better illustrate the method, the MOESP algorithm have to be introduced first to lay the foundation so that the PM-EIV method can be derived. Before that, several assumptions have to be made.Assumption 1.
The difference between the nonlinear functionsf
(
·
), g
(
·
) and their linear approximations respectively, can be neglected.Assumption 2.
Persistent exciting condition is satisfied according to the definition from Ljung [23].Assumption 3.
The system noise sequencesd
(
k
), υ
(
k
), and ω
(
k
) are white noises independent from each other.


### 3.1. MOESP Identification Method

The fundamental feature of MOESP method is to estimate the extended observer matrix of the system based on input/output data. Then system matrices can be obtained through applying least-squares algorithms. To give a brief introduction of MOESP method, system noises are assumed to be zero in this section, that isd
(
k
)
≡
0, υ
(
k
)
≡
0 and ω
(
k
)
≡
0. According to Verhaegen [24], The following equations need to be formed firstly. 								(9)
Y
j
,
i
,
N
=
Γ
i
X
j
,
N
+
H
i
E
j
,
i
,
N
with definitions of Y
j
,
i
,
N, X
j
,
N, and E
j
,
i
,
N as follows:
(10)
Y
j
,
i
,
N
=
[
y
j
⋯
y
j
+
N
-
1
⋮
⋱
⋮
y
j
+
i
-
1
⋯
y
j
+
N
+
i
-
2
]
,
E
j
,
i
,
N
=
[
e
j
⋯
e
j
+
N
-
1
⋮
⋱
⋮
e
j
+
i
-
1
⋯
e
j
+
N
+
i
-
2
]
,
X
j
,
N
=
[
x
j
x
j
+
1
⋯
x
j
+
N
-
1
]
.Therefore, extended observer matrixΓ
i and low triangular block Toeplitz matrix H
i have the following structures:
(11)
Γ
i
=
[
C
T
(
C
A
)
T
⋯
(
C
A
i
-
1
)
T
]
T
,
H
i
=
[
N
0
⋯
0
C
M
N
⋯
0
⋮
⋮
⋱
⋮
C
A
i
-
2
M
C
A
i
-
3
M
⋯
N
]
.Input/output data matrix is formed and factorized by RQ factorization. The following equation can be acquired:(12)
[
E
j
,
i
,
N
Y
j
,
i
,
N
]
=
[
R
11
0
R
21
R
22
]
[
Q
1
Q
2
]
.Combining (12) with (9), we can obtain that:
(13)
Γ
i
X
j
,
N
=
[
R
21
-
H
i
R
11
R
22
]
[
Q
1
Q
2
]
.It can be derived that the column space ofΓ
i and the column space of R
22 are equal. After carrying out SVD of R
22 as in (14), columns of U
n can be regarded as a basis for Γ
i:
(14)
R
22
=
[
U
n
U
n
⊥
]
[
S
1
0
o
S
2
]
[
V
n
T
(
V
n
⊥
)
T
]
.Based on the definition of extended observer matrixΓ
i, transformed system matrices A
T, C
T can be easily calculated with U
n:
(15)
U
n
(
1
)
A
T
=
U
n
(
2
)
,
C
T
=
U
n
(
1
:
l
,
:
)
.AndM
T, N can be acquired by least-squares solutions for. 								(16)
(
U
n
⊥
)
T
H
i
-
(
U
n
⊥
)
T
R
21
R
11
-
1
=
0
.


### 3.2. Errors-in-Variables Problem

Because the identification problem of the Hammerstein AUV model is an EIV one, a more practical subspace identification method handling EIV problem is introduced in this section which takes disturbancesd
(
k
), υ
(
k
), ω
(
k
) into consideration based on the basic MOESP algorithm above. In this case, (9) needs to be modified as:
(17)
Y
~
j
,
i
,
N
=
Γ
i
X
j
,
N
+
H
i
E
~
j
,
i
,
N
-
H
i
V
j
,
i
,
N
′
+
G
i
P
j
,
i
,
N
+
W
j
,
i
,
N
,
where P
j
,
i
,
N, W
j
,
i
,
N are block Hankel matrices of noises d
(
k
) and ω
(
k
). G
i is defined as:
(18)
G
i
=
[
0
0
⋯
0
C
0
⋯
0
⋮
⋮
⋱
⋮
C
A
i
-
2
C
A
i
-
3
⋯
0
]
.Notice thatV
j
,
i
,
N
′ is not only made of υ
(
k
), output measurement noise ω
(
k
) is also involved due to the linearization of nonlinear feedback part. This brings the closed-loop problem into the identification process. One significant obstacle resulting from closed-loop situation is the interference between input signals and output measurement noise which violates the following equation. 								(19)
E
[
u
k
w
j
T
]
=
0
for
(
j
<
k
)
.It means that input signals are no longer unrelated to the past noises. However, problems studied in this paper only need the condition that future noise is unrelated to past input. That is:(20)
E
[
u
k
w
j
T
]
=
0
for
(
j
≥
k
)
.So, in order to eliminate the influence from the noises, instrumental variables can still be formed as past input and output signals. Then following relation can be obtained:(21)
Y
~
i
+
1
,
i
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
=
Γ
i
X
i
+
1
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
+
H
i
E
~
i
+
1
,
i
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
,
(22)
[
E
~
i
+
1
,
i
,
N
E
~
1
,
i
,
N
T
E
~
i
+
1
,
i
,
N
Y
~
1
,
i
,
N
T
Y
~
i
+
1
,
i
,
N
E
~
1
,
i
,
N
T
Y
~
i
+
1
,
i
,
N
Y
~
1
,
i
,
N
T
]
=
[
R
11
(
t
)
0
R
21
(
t
)
R
22
(
t
)
]
[
Q
1
(
t
)
Q
2
(
t
)
]
.According to Theorem  3 of Chou and Verhaegen [13], the column space of Γ
i can be consistently estimated from R
22. A
T, M
T, C
T, and N can be acquired based on OE_PIV algorithm. To save the space, more details can be found in [13].


### 3.3. PM-EIV Subspace Identification

In the above sections, subspace identification method for AUV Hammerstein system has been developed under general noise assumption. As mentioned in Section1, recursive method for AUV identification can be much more suitable due to the properties of oceanic environment and tasks. So in this section, identification method will be modified recursively. One of the key problems handled in recursive identification is the recursive update of SVD process in order to reduce the computational burden which comes with increasing input/output data. In [20], propagator method used in array signal processing area is introduced for recursive subspace identification under PI/PO MOESP schemes. Considering about the EIV problem in the identification of AUV model, propagator method will be extended into EIV framework in this paper for the first time and the resulting algorithm will be called PM-EIV subspace identification method. An important step in PM subspace method is the calculation of observer vector, defined as. 								(23)
z
i
(
t
+
1
)
=
y
~
i
(
t
+
1
)
-
H
i
(
t
+
1
)
e
~
i
(
t
+
1
)
,
where t
=
j
+
N
-
1;  y
~
i
(
t
+
1
), e
~
i
(
t
+
1
) are new updated output/input data vectors respectively.


#### 3.3.1. Update of Observer Vectorz
i in EIV Framework

In this section, RQ factorization method is adopted for updating of observe vector. Based on (22) in Section 3.2, update of the data matrix can be expressed as follows:
(24)
[
E
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
E
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
]
=
[
R
11
(
t
+
1
)
0
R
21
(
t
+
1
)
R
22
(
t
+
1
)
]
[
Q
1
(
t
+
1
)
Q
2
(
t
+
1
)
]
,
where E
~
i
+
1
,
i
,
N
+
1
=
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
],  E
~
1
,
i
,
N
+
1, and Y
~
i
+
1
,
i
,
N
+
1, Y
~
1
,
i
,
N
+
1 are formed in similar manner. A transformation of (24) is as follows:
(25)
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
Y
~
i
+
1
,
i
,
N
y
~
i
(
t
+
1
)
]
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
=
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
[
Q
1
(
t
)
Q
2
(
t
)
I
]
,
where ϕ
(
N
+
1
)
=
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
] denotes the update instrumental variables.Then,given rotation can be implemented to eliminate the e
~
i
(
t
+
1
)
ϕ
(
N
+
1
) in the above equation in order to make a lower triangle form.
(26)
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
Giv
(
t
+
1
)
=
[
R
11
(
t
+
1
)
0
0
R
21
(
t
+
1
)
R
22
(
t
)
z
^
i
(
t
+
1
)
]
.So, the following relation can be obtained:(27)
R
22
(
t
+
1
)
R
22
(
t
+
1
)
T
=
R
22
(
t
)
R
22
(
t
)
T
+
z
^
i
(
t
+
1
)
z
^
i
T
(
t
+
1
)
.In addition, the following equation holds:(28)
R
22
(
t
+
1
)
Q
2
(
t
+
1
)
=
(
Y
~
i
+
1
,
i
,
N
+
1
-
H
^
i
E
~
i
+
1
,
i
,
N
+
1
)
[
E
~
1
,
i
,
N
+
1
T
Y
~
1
,
i
,
N
+
1
T
]
=
R
22
(
t
)
Q
2
(
t
)
+
z
i
(
t
+
1
)
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
.Combining (27) and (28), relation between z
i
(
t
+
1
) and z
^
i
(
t
+
1
) can be found:
(29)
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
=
R
22
(
t
)
Q
2
(
t
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
(
R
22
(
t
)
Q
2
(
t
)
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
.Assuming there is a matrixK
(
t
+
1
) satisfies the following equation:
(30)
K
(
t
+
1
)
K
(
t
+
1
)
T
=
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
+
R
22
(
t
)
R
22
(
t
)
T
.Then observer vectorz
i can be updated according to (31):
(31)
z
i
(
t
+
1
)
[
e
i
T
(
N
+
1
)
y
i
T
(
N
+
1
)
]
=
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
,
z
i
(
t
+
1
)
=
K
ϕ
-
1
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
ϕ
N
T
,
where K
ϕ
=
ϕ
N
ϕ
N
T is the coefficient related to instrumental variables, and notice that K
(
t
+
1
) is not necessary to be square.


#### 3.3.2. Estimation of Observer MatrixΓ
i

Since the observer vectorz
i can be obtained in EIV scheme. Then, extended observer matrix Γ
i can be found through propagator method. Observer matrix Γ
i can be expressed in the following form. 										(32)
Γ
i
=
[
Γ
i
1
Γ
i
2
]
=
[
I
P
m
T
]
Γ
i
1
=
Q
s
Γ
i
1
.Since the matrixΓ
i
1
∈
R
n
×
n has full rank, column space of Γ
i equals that of Q
s. Combining with (23), observer vector z
i can be divided as:
(33)
z
i
(
t
+
1
)
=
Q
s
Γ
i
1
x
(
t
+
1
)
+
b
i
(
t
+
1
)
=
[
z
i
1
(
t
+
1
)
z
i
2
(
t
+
1
)
]
.Without consideration of noise termb
i
(
t
+
1
), it can be easily established that z
i
2
=
P
m
T
z
i
1. Then the P
m
T can be solved with least-squares methods. However, existence of noise term will lead to a biased estimation of P
m
T. So the IVPM algorithm proposed by Mercere [20] is adopted to estimate P
m
T. A suitable variable γ
∈
R
n
×
1 needs to be found with no correlation with system noises. According to Section 3.2, past system input date may satisfy the condition.Finally, recursive estimation ofP
m
T can be acquired in the following RLS form:
(34)
K
p
(
t
+
1
)
=
γ
T
(
t
+
1
)
R
(
t
)
(
1
+
γ
T
(
t
+
1
)
R
(
t
)
z
i
1
(
t
+
1
)
)
-
1
,
P
m
T
(
t
+
1
)
=
P
m
T
(
t
)
+
[
z
i
2
(
t
+
1
)
-
P
m
T
(
t
)
z
i
1
(
t
+
1
)
]
K
p
(
t
+
1
)
,
R
(
t
+
1
)
=
R
(
t
)
-
R
(
t
)
z
i
1
(
t
+
1
)
K
p
(
t
+
1
)
,
where R
(
t
)
=
E
[
z
i
1
(
t
)
γ
T
(
t
)
]
-
1.With the estimation of observer matrixΓ
i, algorithms mentioned in Section 3.2 can be applied here to estimate the system matrices A
T, M
T, C
T, and N. Therefore, propagator based subspace identification method in errors-in-variables scheme (PM-EIV) is derived. Based on the IVPM method, the algorithm is more suitable for the identification of Hammerstein AUV model proposed in Section 2. In the following, the identification algorithm will be verified through experiments carried out by the AUV simulation platform.


## 3.1. MOESP Identification Method

The fundamental feature of MOESP method is to estimate the extended observer matrix of the system based on input/output data. Then system matrices can be obtained through applying least-squares algorithms. To give a brief introduction of MOESP method, system noises are assumed to be zero in this section, that isd
(
k
)
≡
0, υ
(
k
)
≡
0 and ω
(
k
)
≡
0. According to Verhaegen [24], The following equations need to be formed firstly. 								(9)
Y
j
,
i
,
N
=
Γ
i
X
j
,
N
+
H
i
E
j
,
i
,
N
with definitions of Y
j
,
i
,
N, X
j
,
N, and E
j
,
i
,
N as follows:
(10)
Y
j
,
i
,
N
=
[
y
j
⋯
y
j
+
N
-
1
⋮
⋱
⋮
y
j
+
i
-
1
⋯
y
j
+
N
+
i
-
2
]
,
E
j
,
i
,
N
=
[
e
j
⋯
e
j
+
N
-
1
⋮
⋱
⋮
e
j
+
i
-
1
⋯
e
j
+
N
+
i
-
2
]
,
X
j
,
N
=
[
x
j
x
j
+
1
⋯
x
j
+
N
-
1
]
.Therefore, extended observer matrixΓ
i and low triangular block Toeplitz matrix H
i have the following structures:
(11)
Γ
i
=
[
C
T
(
C
A
)
T
⋯
(
C
A
i
-
1
)
T
]
T
,
H
i
=
[
N
0
⋯
0
C
M
N
⋯
0
⋮
⋮
⋱
⋮
C
A
i
-
2
M
C
A
i
-
3
M
⋯
N
]
.Input/output data matrix is formed and factorized by RQ factorization. The following equation can be acquired:(12)
[
E
j
,
i
,
N
Y
j
,
i
,
N
]
=
[
R
11
0
R
21
R
22
]
[
Q
1
Q
2
]
.Combining (12) with (9), we can obtain that:
(13)
Γ
i
X
j
,
N
=
[
R
21
-
H
i
R
11
R
22
]
[
Q
1
Q
2
]
.It can be derived that the column space ofΓ
i and the column space of R
22 are equal. After carrying out SVD of R
22 as in (14), columns of U
n can be regarded as a basis for Γ
i:
(14)
R
22
=
[
U
n
U
n
⊥
]
[
S
1
0
o
S
2
]
[
V
n
T
(
V
n
⊥
)
T
]
.Based on the definition of extended observer matrixΓ
i, transformed system matrices A
T, C
T can be easily calculated with U
n:
(15)
U
n
(
1
)
A
T
=
U
n
(
2
)
,
C
T
=
U
n
(
1
:
l
,
:
)
.AndM
T, N can be acquired by least-squares solutions for. 								(16)
(
U
n
⊥
)
T
H
i
-
(
U
n
⊥
)
T
R
21
R
11
-
1
=
0
.


## 3.2. Errors-in-Variables Problem

Because the identification problem of the Hammerstein AUV model is an EIV one, a more practical subspace identification method handling EIV problem is introduced in this section which takes disturbancesd
(
k
), υ
(
k
), ω
(
k
) into consideration based on the basic MOESP algorithm above. In this case, (9) needs to be modified as:
(17)
Y
~
j
,
i
,
N
=
Γ
i
X
j
,
N
+
H
i
E
~
j
,
i
,
N
-
H
i
V
j
,
i
,
N
′
+
G
i
P
j
,
i
,
N
+
W
j
,
i
,
N
,
where P
j
,
i
,
N, W
j
,
i
,
N are block Hankel matrices of noises d
(
k
) and ω
(
k
). G
i is defined as:
(18)
G
i
=
[
0
0
⋯
0
C
0
⋯
0
⋮
⋮
⋱
⋮
C
A
i
-
2
C
A
i
-
3
⋯
0
]
.Notice thatV
j
,
i
,
N
′ is not only made of υ
(
k
), output measurement noise ω
(
k
) is also involved due to the linearization of nonlinear feedback part. This brings the closed-loop problem into the identification process. One significant obstacle resulting from closed-loop situation is the interference between input signals and output measurement noise which violates the following equation. 								(19)
E
[
u
k
w
j
T
]
=
0
for
(
j
<
k
)
.It means that input signals are no longer unrelated to the past noises. However, problems studied in this paper only need the condition that future noise is unrelated to past input. That is:(20)
E
[
u
k
w
j
T
]
=
0
for
(
j
≥
k
)
.So, in order to eliminate the influence from the noises, instrumental variables can still be formed as past input and output signals. Then following relation can be obtained:(21)
Y
~
i
+
1
,
i
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
=
Γ
i
X
i
+
1
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
+
H
i
E
~
i
+
1
,
i
,
N
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
]
,
(22)
[
E
~
i
+
1
,
i
,
N
E
~
1
,
i
,
N
T
E
~
i
+
1
,
i
,
N
Y
~
1
,
i
,
N
T
Y
~
i
+
1
,
i
,
N
E
~
1
,
i
,
N
T
Y
~
i
+
1
,
i
,
N
Y
~
1
,
i
,
N
T
]
=
[
R
11
(
t
)
0
R
21
(
t
)
R
22
(
t
)
]
[
Q
1
(
t
)
Q
2
(
t
)
]
.According to Theorem  3 of Chou and Verhaegen [13], the column space of Γ
i can be consistently estimated from R
22. A
T, M
T, C
T, and N can be acquired based on OE_PIV algorithm. To save the space, more details can be found in [13].


## 3.3. PM-EIV Subspace Identification

In the above sections, subspace identification method for AUV Hammerstein system has been developed under general noise assumption. As mentioned in Section1, recursive method for AUV identification can be much more suitable due to the properties of oceanic environment and tasks. So in this section, identification method will be modified recursively. One of the key problems handled in recursive identification is the recursive update of SVD process in order to reduce the computational burden which comes with increasing input/output data. In [20], propagator method used in array signal processing area is introduced for recursive subspace identification under PI/PO MOESP schemes. Considering about the EIV problem in the identification of AUV model, propagator method will be extended into EIV framework in this paper for the first time and the resulting algorithm will be called PM-EIV subspace identification method. An important step in PM subspace method is the calculation of observer vector, defined as. 								(23)
z
i
(
t
+
1
)
=
y
~
i
(
t
+
1
)
-
H
i
(
t
+
1
)
e
~
i
(
t
+
1
)
,
where t
=
j
+
N
-
1;  y
~
i
(
t
+
1
), e
~
i
(
t
+
1
) are new updated output/input data vectors respectively.


### 3.3.1. Update of Observer Vectorz
i in EIV Framework

In this section, RQ factorization method is adopted for updating of observe vector. Based on (22) in Section 3.2, update of the data matrix can be expressed as follows:
(24)
[
E
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
E
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
]
=
[
R
11
(
t
+
1
)
0
R
21
(
t
+
1
)
R
22
(
t
+
1
)
]
[
Q
1
(
t
+
1
)
Q
2
(
t
+
1
)
]
,
where E
~
i
+
1
,
i
,
N
+
1
=
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
],  E
~
1
,
i
,
N
+
1, and Y
~
i
+
1
,
i
,
N
+
1, Y
~
1
,
i
,
N
+
1 are formed in similar manner. A transformation of (24) is as follows:
(25)
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
Y
~
i
+
1
,
i
,
N
y
~
i
(
t
+
1
)
]
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
=
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
[
Q
1
(
t
)
Q
2
(
t
)
I
]
,
where ϕ
(
N
+
1
)
=
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
] denotes the update instrumental variables.Then,given rotation can be implemented to eliminate the e
~
i
(
t
+
1
)
ϕ
(
N
+
1
) in the above equation in order to make a lower triangle form.
(26)
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
Giv
(
t
+
1
)
=
[
R
11
(
t
+
1
)
0
0
R
21
(
t
+
1
)
R
22
(
t
)
z
^
i
(
t
+
1
)
]
.So, the following relation can be obtained:(27)
R
22
(
t
+
1
)
R
22
(
t
+
1
)
T
=
R
22
(
t
)
R
22
(
t
)
T
+
z
^
i
(
t
+
1
)
z
^
i
T
(
t
+
1
)
.In addition, the following equation holds:(28)
R
22
(
t
+
1
)
Q
2
(
t
+
1
)
=
(
Y
~
i
+
1
,
i
,
N
+
1
-
H
^
i
E
~
i
+
1
,
i
,
N
+
1
)
[
E
~
1
,
i
,
N
+
1
T
Y
~
1
,
i
,
N
+
1
T
]
=
R
22
(
t
)
Q
2
(
t
)
+
z
i
(
t
+
1
)
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
.Combining (27) and (28), relation between z
i
(
t
+
1
) and z
^
i
(
t
+
1
) can be found:
(29)
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
=
R
22
(
t
)
Q
2
(
t
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
(
R
22
(
t
)
Q
2
(
t
)
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
.Assuming there is a matrixK
(
t
+
1
) satisfies the following equation:
(30)
K
(
t
+
1
)
K
(
t
+
1
)
T
=
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
+
R
22
(
t
)
R
22
(
t
)
T
.Then observer vectorz
i can be updated according to (31):
(31)
z
i
(
t
+
1
)
[
e
i
T
(
N
+
1
)
y
i
T
(
N
+
1
)
]
=
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
,
z
i
(
t
+
1
)
=
K
ϕ
-
1
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
ϕ
N
T
,
where K
ϕ
=
ϕ
N
ϕ
N
T is the coefficient related to instrumental variables, and notice that K
(
t
+
1
) is not necessary to be square.


### 3.3.2. Estimation of Observer MatrixΓ
i

Since the observer vectorz
i can be obtained in EIV scheme. Then, extended observer matrix Γ
i can be found through propagator method. Observer matrix Γ
i can be expressed in the following form. 										(32)
Γ
i
=
[
Γ
i
1
Γ
i
2
]
=
[
I
P
m
T
]
Γ
i
1
=
Q
s
Γ
i
1
.Since the matrixΓ
i
1
∈
R
n
×
n has full rank, column space of Γ
i equals that of Q
s. Combining with (23), observer vector z
i can be divided as:
(33)
z
i
(
t
+
1
)
=
Q
s
Γ
i
1
x
(
t
+
1
)
+
b
i
(
t
+
1
)
=
[
z
i
1
(
t
+
1
)
z
i
2
(
t
+
1
)
]
.Without consideration of noise termb
i
(
t
+
1
), it can be easily established that z
i
2
=
P
m
T
z
i
1. Then the P
m
T can be solved with least-squares methods. However, existence of noise term will lead to a biased estimation of P
m
T. So the IVPM algorithm proposed by Mercere [20] is adopted to estimate P
m
T. A suitable variable γ
∈
R
n
×
1 needs to be found with no correlation with system noises. According to Section 3.2, past system input date may satisfy the condition.Finally, recursive estimation ofP
m
T can be acquired in the following RLS form:
(34)
K
p
(
t
+
1
)
=
γ
T
(
t
+
1
)
R
(
t
)
(
1
+
γ
T
(
t
+
1
)
R
(
t
)
z
i
1
(
t
+
1
)
)
-
1
,
P
m
T
(
t
+
1
)
=
P
m
T
(
t
)
+
[
z
i
2
(
t
+
1
)
-
P
m
T
(
t
)
z
i
1
(
t
+
1
)
]
K
p
(
t
+
1
)
,
R
(
t
+
1
)
=
R
(
t
)
-
R
(
t
)
z
i
1
(
t
+
1
)
K
p
(
t
+
1
)
,
where R
(
t
)
=
E
[
z
i
1
(
t
)
γ
T
(
t
)
]
-
1.With the estimation of observer matrixΓ
i, algorithms mentioned in Section 3.2 can be applied here to estimate the system matrices A
T, M
T, C
T, and N. Therefore, propagator based subspace identification method in errors-in-variables scheme (PM-EIV) is derived. Based on the IVPM method, the algorithm is more suitable for the identification of Hammerstein AUV model proposed in Section 2. In the following, the identification algorithm will be verified through experiments carried out by the AUV simulation platform.


## 3.3.1. Update of Observer Vectorz
i in EIV Framework

In this section, RQ factorization method is adopted for updating of observe vector. Based on (22) in Section 3.2, update of the data matrix can be expressed as follows:
(24)
[
E
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
E
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
E
~
1
,
i
,
N
+
1
T
Y
~
i
+
1
,
i
,
N
+
1
Y
~
1
,
i
,
N
+
1
T
]
=
[
R
11
(
t
+
1
)
0
R
21
(
t
+
1
)
R
22
(
t
+
1
)
]
[
Q
1
(
t
+
1
)
Q
2
(
t
+
1
)
]
,
where E
~
i
+
1
,
i
,
N
+
1
=
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
],  E
~
1
,
i
,
N
+
1, and Y
~
i
+
1
,
i
,
N
+
1, Y
~
1
,
i
,
N
+
1 are formed in similar manner. A transformation of (24) is as follows:
(25)
[
E
~
i
+
1
,
i
,
N
e
~
i
(
t
+
1
)
Y
~
i
+
1
,
i
,
N
y
~
i
(
t
+
1
)
]
[
E
~
1
,
i
,
N
T
Y
~
1
,
i
,
N
T
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
=
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
[
Q
1
(
t
)
Q
2
(
t
)
I
]
,
where ϕ
(
N
+
1
)
=
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
] denotes the update instrumental variables.Then,given rotation can be implemented to eliminate the e
~
i
(
t
+
1
)
ϕ
(
N
+
1
) in the above equation in order to make a lower triangle form.
(26)
[
R
11
(
t
)
0
e
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
R
21
(
t
)
R
22
(
t
)
y
~
i
(
t
+
1
)
ϕ
(
N
+
1
)
]
Giv
(
t
+
1
)
=
[
R
11
(
t
+
1
)
0
0
R
21
(
t
+
1
)
R
22
(
t
)
z
^
i
(
t
+
1
)
]
.So, the following relation can be obtained:(27)
R
22
(
t
+
1
)
R
22
(
t
+
1
)
T
=
R
22
(
t
)
R
22
(
t
)
T
+
z
^
i
(
t
+
1
)
z
^
i
T
(
t
+
1
)
.In addition, the following equation holds:(28)
R
22
(
t
+
1
)
Q
2
(
t
+
1
)
=
(
Y
~
i
+
1
,
i
,
N
+
1
-
H
^
i
E
~
i
+
1
,
i
,
N
+
1
)
[
E
~
1
,
i
,
N
+
1
T
Y
~
1
,
i
,
N
+
1
T
]
=
R
22
(
t
)
Q
2
(
t
)
+
z
i
(
t
+
1
)
[
e
~
i
(
N
+
1
)
T
y
~
i
(
N
+
1
)
T
]
.Combining (27) and (28), relation between z
i
(
t
+
1
) and z
^
i
(
t
+
1
) can be found:
(29)
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
=
R
22
(
t
)
Q
2
(
t
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
(
R
22
(
t
)
Q
2
(
t
)
)
T
+
z
i
(
t
+
1
)
ϕ
(
N
+
1
)
ϕ
(
N
+
1
)
T
z
i
(
t
+
1
)
T
.Assuming there is a matrixK
(
t
+
1
) satisfies the following equation:
(30)
K
(
t
+
1
)
K
(
t
+
1
)
T
=
z
^
i
(
t
+
1
)
z
^
i
(
t
+
1
)
T
+
R
22
(
t
)
R
22
(
t
)
T
.Then observer vectorz
i can be updated according to (31):
(31)
z
i
(
t
+
1
)
[
e
i
T
(
N
+
1
)
y
i
T
(
N
+
1
)
]
=
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
,
z
i
(
t
+
1
)
=
K
ϕ
-
1
(
K
(
t
+
1
)
-
R
22
(
t
)
Q
2
(
t
)
)
ϕ
N
T
,
where K
ϕ
=
ϕ
N
ϕ
N
T is the coefficient related to instrumental variables, and notice that K
(
t
+
1
) is not necessary to be square.


## 3.3.2. Estimation of Observer MatrixΓ
i

Since the observer vectorz
i can be obtained in EIV scheme. Then, extended observer matrix Γ
i can be found through propagator method. Observer matrix Γ
i can be expressed in the following form. 										(32)
Γ
i
=
[
Γ
i
1
Γ
i
2
]
=
[
I
P
m
T
]
Γ
i
1
=
Q
s
Γ
i
1
.Since the matrixΓ
i
1
∈
R
n
×
n has full rank, column space of Γ
i equals that of Q
s. Combining with (23), observer vector z
i can be divided as:
(33)
z
i
(
t
+
1
)
=
Q
s
Γ
i
1
x
(
t
+
1
)
+
b
i
(
t
+
1
)
=
[
z
i
1
(
t
+
1
)
z
i
2
(
t
+
1
)
]
.Without consideration of noise termb
i
(
t
+
1
), it can be easily established that z
i
2
=
P
m
T
z
i
1. Then the P
m
T can be solved with least-squares methods. However, existence of noise term will lead to a biased estimation of P
m
T. So the IVPM algorithm proposed by Mercere [20] is adopted to estimate P
m
T. A suitable variable γ
∈
R
n
×
1 needs to be found with no correlation with system noises. According to Section 3.2, past system input date may satisfy the condition.Finally, recursive estimation ofP
m
T can be acquired in the following RLS form:
(34)
K
p
(
t
+
1
)
=
γ
T
(
t
+
1
)
R
(
t
)
(
1
+
γ
T
(
t
+
1
)
R
(
t
)
z
i
1
(
t
+
1
)
)
-
1
,
P
m
T
(
t
+
1
)
=
P
m
T
(
t
)
+
[
z
i
2
(
t
+
1
)
-
P
m
T
(
t
)
z
i
1
(
t
+
1
)
]
K
p
(
t
+
1
)
,
R
(
t
+
1
)
=
R
(
t
)
-
R
(
t
)
z
i
1
(
t
+
1
)
K
p
(
t
+
1
)
,
where R
(
t
)
=
E
[
z
i
1
(
t
)
γ
T
(
t
)
]
-
1.With the estimation of observer matrixΓ
i, algorithms mentioned in Section 3.2 can be applied here to estimate the system matrices A
T, M
T, C
T, and N. Therefore, propagator based subspace identification method in errors-in-variables scheme (PM-EIV) is derived. Based on the IVPM method, the algorithm is more suitable for the identification of Hammerstein AUV model proposed in Section 2. In the following, the identification algorithm will be verified through experiments carried out by the AUV simulation platform.


## 4. Simulations and Results

In this section, simulation experiments are carried out to evaluate the performance of the PM-EIV algorithm proposed in this paper through identifying the Hammerstein AUV model based on data from the AUV simulation platform shown as in Figures2-3. After a brief introduction of the AUV simulation platform, two typical identification cases are investigated. One is the identification of AUV model based on the MOESP method without consideration of noise. The other one is the verification of PM-EIV algorithm under general noise assumption. Finally, to be more practical, model identification based on data from a closed-loop path following simulation experiment is performed.Figure 2
Information flow of AUV simulation platform.Figure 3
AUV simulation platform.


### 4.1. AUV Simulation Platform

The basic structure of AUV simulation platform is depicted as in Figure3. The whole system is connected through Ethernet and responsibilities of five main components are introduced here.(C1) Surface Interface Computer (SIC). Surface interface software is running on SIC which is in charge of deploying missions for AUV and monitoring the states of the system. The software also allows for manual intervention in case of emergency.(C2) Mission Control Computer (MiC). Mission Management software developed in QNX real-time system is the core of MiC. MiC is mainly responsible for path plan according to missions from SIC and fault diagnosis of the system.(C3) Motion Control Computer (MoC). Motion control software running on MoC aims at controlling the states such as heading, speed and depth of AUV based on the preplanned paths from MiC. In field experiments, MoC is also in charge of navigation of the system.(C4) Model Computer (MC). MC is the host for mathematic models of AUV, thruster and control surface. The mathematic AUV model is a validated model with full coefficients obtained from water-tank experiments.(C5) 3D Simulation Computer (3DC). 3DC provides an approach for visual simulation relied on vir-tual reality technology. Vega Prime is used to construct the virtual oceanic environment and AUV model is developed in Multigen Creator.A more detailed process flow is described in Figure2.


### 4.2. Case  1: Identification without Consideration of Noise

In an ideal situation, process noise and measurement noises can be ignored, so that the ordinary MOESP algorithm can be applied to identify the Hammerstein AUV model described as in Figure1. This case aims at testifying the reasonability of the Hammerstein model for an AUV system. All system input signals are chosen as sinusoid curves with different amplitudes and periods and the value of N in (4) is set to be 4. System identification results based on O-MOESP algorithm can be obtained. Figures 4, 5, and 6 shows the prediction errors between the outputs of the identified model and the original outputs. Three main system outputs surge velocity, pitch rate, and yaw rate, that play important roles in control and navigation of underactuated AUV which are considered here. From Figures 4–6, it can be concluded that even though identification errors exist, Hammerstein model constructed in Section 2 can act as a suitable structure for AUV dynamics.Figure 4
Prediction of surge speed.Figure 5
Prediction of yaw rate.Figure 6
Prediction of pitch rate.


### 4.3. Case  2: Identification under EIV Framework

In this case, general noises are added on the model computer in the AUV simulation platform. So the identification problem becomes an EIV one which can be solved recursively by the PM-EIV algorithm proposed in this paper. System inputs of the simulation platform are still sinusoid curves and the value ofN is chosen to be 6. Since the PM-EIV algorithm is a recursive one, a relative small amount of data is used to acquire an initial value of the model at the beginning. Identification results are shown in Figures 7, 8, and 9. Based on the Figures 7–9, it is reliable to conclude that the PM-EIV algorithm is effective and feasible in identifying the Hammerstein AUV model under general noise assumption.Figure 7
Prediction of surge speed.Figure 8
Prediction of system pitch rate.Figure 9
Prediction of system yaw rate.Remark. Through the above two different identification cases under different situations, it has been approved that Hammerstein model proposed in Section 2 is capable of representing the AUV dynamic system and the PM-EIV method is able to identify the Hammerstein model recursively under general noise assumption. It is also interesting to notice that steady state prediction errors in Case 2 are decreased comparison with those in Case 1, it is also interesting to notice that steady state prediction errors in Case 2 are decreased comparing with those in Case 1 due to the increase of N adopted.


### 4.4. Identification from Closed-Loop Simulation

In the above simulations, AUV model is identified using data from open loop control of surge speed, heading and yaw. However, in practical applications, data for model identification is usually collected from field experiments with close loop control for specific preplanned paths. So in this section, AUV identification process is carried out based on data from a closed-loop simulation experiment. To be more pellucid, the feasibility of PM-EIV algorithm is not illustrated by predicting the surge speed, pitch rate and yaw rate, but shown by predicting the trajectory of AUV in the test. Figures10 and 11 are the simulation results from the AUV simulation platform. The preplanned path is a circle with an origin at (500, 500) m and radius of 300 m. The depth command is 10 m. Simulation results have shown that AUV can follow the preplanned circle very well and the diving process is stable and fast.Simulation results: (a) horizontal projection of the trajectory; (b) 3D trajectory of AUV.
(a)
(b)Figure 11
Identification result.Then identification of the Hammerstein AUV model is based on the inputs/outputs of this experiment. Noises are also considered. Figure11 has shown the identification result in a different point of view. It can be seen that prediction error between the identified model of AUV and the actual trajectory of AUV is large at the beginning because the recursive identification procedure needs time to converge. Therefore, a position calibration operation is carried out when the identification results have converged at simulation time 350 s. The Cali-Point in Figure 11 is where the calibration is implemented. The green line indicates the distance between Cali-Point and expected point. After the calibration, the trajectory of the identified Hammerstein AUV model can follow the trajectory of AUV with satisfying accuracy.


## 4.1. AUV Simulation Platform

The basic structure of AUV simulation platform is depicted as in Figure3. The whole system is connected through Ethernet and responsibilities of five main components are introduced here.(C1) Surface Interface Computer (SIC). Surface interface software is running on SIC which is in charge of deploying missions for AUV and monitoring the states of the system. The software also allows for manual intervention in case of emergency.(C2) Mission Control Computer (MiC). Mission Management software developed in QNX real-time system is the core of MiC. MiC is mainly responsible for path plan according to missions from SIC and fault diagnosis of the system.(C3) Motion Control Computer (MoC). Motion control software running on MoC aims at controlling the states such as heading, speed and depth of AUV based on the preplanned paths from MiC. In field experiments, MoC is also in charge of navigation of the system.(C4) Model Computer (MC). MC is the host for mathematic models of AUV, thruster and control surface. The mathematic AUV model is a validated model with full coefficients obtained from water-tank experiments.(C5) 3D Simulation Computer (3DC). 3DC provides an approach for visual simulation relied on vir-tual reality technology. Vega Prime is used to construct the virtual oceanic environment and AUV model is developed in Multigen Creator.A more detailed process flow is described in Figure2.


## 4.2. Case  1: Identification without Consideration of Noise

In an ideal situation, process noise and measurement noises can be ignored, so that the ordinary MOESP algorithm can be applied to identify the Hammerstein AUV model described as in Figure1. This case aims at testifying the reasonability of the Hammerstein model for an AUV system. All system input signals are chosen as sinusoid curves with different amplitudes and periods and the value of N in (4) is set to be 4. System identification results based on O-MOESP algorithm can be obtained. Figures 4, 5, and 6 shows the prediction errors between the outputs of the identified model and the original outputs. Three main system outputs surge velocity, pitch rate, and yaw rate, that play important roles in control and navigation of underactuated AUV which are considered here. From Figures 4–6, it can be concluded that even though identification errors exist, Hammerstein model constructed in Section 2 can act as a suitable structure for AUV dynamics.Figure 4
Prediction of surge speed.Figure 5
Prediction of yaw rate.Figure 6
Prediction of pitch rate.


## 4.3. Case  2: Identification under EIV Framework

In this case, general noises are added on the model computer in the AUV simulation platform. So the identification problem becomes an EIV one which can be solved recursively by the PM-EIV algorithm proposed in this paper. System inputs of the simulation platform are still sinusoid curves and the value ofN is chosen to be 6. Since the PM-EIV algorithm is a recursive one, a relative small amount of data is used to acquire an initial value of the model at the beginning. Identification results are shown in Figures 7, 8, and 9. Based on the Figures 7–9, it is reliable to conclude that the PM-EIV algorithm is effective and feasible in identifying the Hammerstein AUV model under general noise assumption.Figure 7
Prediction of surge speed.Figure 8
Prediction of system pitch rate.Figure 9
Prediction of system yaw rate.Remark. Through the above two different identification cases under different situations, it has been approved that Hammerstein model proposed in Section 2 is capable of representing the AUV dynamic system and the PM-EIV method is able to identify the Hammerstein model recursively under general noise assumption. It is also interesting to notice that steady state prediction errors in Case 2 are decreased comparison with those in Case 1, it is also interesting to notice that steady state prediction errors in Case 2 are decreased comparing with those in Case 1 due to the increase of N adopted.


## 4.4. Identification from Closed-Loop Simulation

In the above simulations, AUV model is identified using data from open loop control of surge speed, heading and yaw. However, in practical applications, data for model identification is usually collected from field experiments with close loop control for specific preplanned paths. So in this section, AUV identification process is carried out based on data from a closed-loop simulation experiment. To be more pellucid, the feasibility of PM-EIV algorithm is not illustrated by predicting the surge speed, pitch rate and yaw rate, but shown by predicting the trajectory of AUV in the test. Figures10 and 11 are the simulation results from the AUV simulation platform. The preplanned path is a circle with an origin at (500, 500) m and radius of 300 m. The depth command is 10 m. Simulation results have shown that AUV can follow the preplanned circle very well and the diving process is stable and fast.Simulation results: (a) horizontal projection of the trajectory; (b) 3D trajectory of AUV.
(a)
(b)Figure 11
Identification result.Then identification of the Hammerstein AUV model is based on the inputs/outputs of this experiment. Noises are also considered. Figure11 has shown the identification result in a different point of view. It can be seen that prediction error between the identified model of AUV and the actual trajectory of AUV is large at the beginning because the recursive identification procedure needs time to converge. Therefore, a position calibration operation is carried out when the identification results have converged at simulation time 350 s. The Cali-Point in Figure 11 is where the calibration is implemented. The green line indicates the distance between Cali-Point and expected point. After the calibration, the trajectory of the identified Hammerstein AUV model can follow the trajectory of AUV with satisfying accuracy.


## 5. Conclusions

In this paper, a recursive subspace identification algorithm PM-EIV is derived under general noise assumption. subspace identification algorithm based on propagator method is extended into EIV framework is extended into EIV framework. In order to implement the method on identification of AUV model with consideration about nonlinearities and couplings at the same time, a Hammerstein AUV model is constructed for the first time. Three simulation experiments under different conditions are carried out to verify the feasibility of the model and the effectiveness of the proposed algorithm.In the future study, system noises will not be restricted to white noise and noise models related with oceanic environment and sensors will be introduced to make the algorithm more practical.



---

*Source: 547539-2014-01-16.xml*