# Application of Recursive Subspace Method in Vehicle Lateral Dynamics Model Identification

Authors: Tengyue Ba; Xiqiang Guan; Jian W. Zhang; Sanzhou Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Modeling of vehicle behavior based on the identification method has received a renewed attention in recent years. In order to improve the linear time-invariant vehicle identification model, a more general identifiable vehicle model structure is proposed, in which time-varying characteristics of vehicle speed and cornering stiffness are taken into consideration. To identify the proposed linear time-varying vehicle model, a well-established data-driven method, named recursive optimized version of predictor-based subspace identification, is introduced. Before vehicle model identification, the influences of four parameters in the subspace algorithm are studied based on pulse input road test. And then a set of practical optimal parameters are chosen and used for the vehicle model identification. Through a series of standard road tests under different maneuvers, the linear time-varying vehicle model can be identified in real-time. It is demonstrated by comparison that the predicted outputs of the proposed vehicle model are much closer to the real vehicle outputs and the model has a wider range of application.


---

## Body

## 1. Introduction

With the development of vehicle system and its control, the modeling of vehicle dynamics behavior has received a renewed attention in recent years. On the one hand, the vehicle model should be as precise as possible to describe the real vehicle performance under different conditions. On the other hand, the model should not be too complicated in consideration of time and cost. Generally, the model derived from physical law is the first choice, because it can provide both physical insight and good control basis. A well-known linear time-invariant (LTI) vehicle model, called single track model [1, 2], is widely used in vehicle system analysis and controller design [3, 4]. As pointed out by Canale et al. [5], the real vehicle response is “too different” from the employed model, which suffers from undermodeling and parameter uncertainties. Fortunately, the modeling method based on identification can provide a good tradeoff between model complexity and accuracy; and it has been widely used in vehicle systems [6].A linear parameter varying (LPV) vehicle model [7, 8] is proposed for vehicle model identification. In the model structure, the coefficients are assumed to be polynomial function of parameters a
y or 1/v
x. To identify the multi-input multioutput (MIMO) vehicle system, Cerone et al. [8] addressed a set-membership method, but the complex and low efficiency identification process consisted of two single-input single-output (SISO) LPV representations. Moreover, the model structure is not practical, because the sideslip angle chosen as the output signal cannot be measured directly for most real vehicles. Guan et al. [9] presented an identifiable LTI vehicle model structure and applied subspace identification methods (SIMs) to estimate the vehicle system. In their model structure, the measurable lateral acceleration was used as the output signal, but the model was limited to a given vehicle speed and linear range for lateral acceleration. Considering the time-varying characteristics of vehicle speed and cornering stiffness, a more general vehicle model is needed. Therefore, a major contribution of this work is to propose a low-order identifiable model structure that describes the real vehicle behavior, for the measurable lateral acceleration serves as an output signal. This can be more practical and the LTV vehicle model structure is more general.Another contribution is the application of practical subspace identification to vehicle modeling. In fact, the traditional SIMs, including CVA, N4SID, and MOESP [10, 11], are successfully used in real system [12], and they do not suffer from the nonlinear optimization problem because of the robust numerical tools, such as the RQ factorization and the singular values decomposition (SVD). But they usually give a biased estimation in closed-loop operation and are not suitable for LTV system. In order to overcome the deficiencies of the traditional SIMs, Chiuso [13] proposed a subspace method named P
B
S
I
D
o
p
t, which can provide asymptotically consistent estimation even in closed-loop operations and is currently one of the most promising solutions [14]. On the other hand, a recursive form of subspace method (RSIM) is necessary for the problem of LTV model structure, and there are two ways to implement it. One way [15] was to improve the computational speed and complexity of SVD, but these technologies were limited in practice, because the disturbance was strictly required as temporally and spatially white noise. The other way was to update the extended observability matrix directly so as to avoid the SVD. Kameyama et al. [16] presented a local RSIM algorithm with fixed data size, but the computation was too heavy. Mercère et al. [17] proposed a propagator method to estimate state vector, but it was complicated to find a special permutation matrix before identification. As mentioned in literature [18], a bad choice of a permutation matrix might result in a wrong performance of identification. Based on Yang’s criterion [19], projection approximation subspace tracking technology (PAST) [20] was used to estimate the state sequence, and the method could give a consistent estimation in spite of colored process noises. As mentioned in literature [21], the PAST may oscillate without convergence, and an improvement of orthonormal PAST (OPAST) was put forward, which had a global convergence property. In this work, a recursive version of P
B
S
I
D
o
p
t method based on OPAST is proposed for the purpose of vehicle model identification.In order to validate the proposed vehicle identification model, a series of vehicle standard road tests are carried out. Before the vehicle model identification, the effect of four parameters in subspace algorithm is studied based on pulse input test, and a set of optimal parameters is determined for the tested vehicle. By comparison with LTI identification model, it will be shown that the improved LTV model can describe vehicle behavior precisely and enlarge the model’s application scope. Moreover, the proposed subspace identification method is proved practical and effective.


## 2. Vehicle Model Identification Method

### 2.1. Identifiable Vehicle Model Structure

A commonly used model in vehicle lateral dynamics research is the single track model, which can simulate the lateral behavior under many conditions. Under the assumptions of small steering angle and same slip angles of left and right tires, and ignoring the roll motion, the differential equations of vehicle lateral motion and yaw motion are expressed as below:(1)
m
v
˙
y
+
v
x
γ
=
F
f
+
F
r
,
I
z
γ
˙
=
a
F
f
-
b
F
r
,where m is the vehicle mass, v
x is the longitudinal speed, v
y is the lateral speed of the CG of the vehicle, γ is the yaw rate, I
z is the vehicle yaw moment inertia, and F
f and F
r are front and rear tire lateral forces (Figure 1).Figure 1
The 2-DOF single track model.For linear tire model under small tire slip angle, the lateral forces can be written as(2)
F
f
=
-
K
f
β
f
,
F
r
=
-
K
r
β
r
,where K
f and K
r are the average cornering stiffness of front and rear axles and β
f and β
r are the front and rear tire slip angles, respectively. The tire slip angles can be expressed approximatively as(3)
β
f
=
-
δ
+
v
y
+
a
γ
v
x
,
β
r
=
v
y
-
b
γ
v
x
.Substituting (2) and (3) into (1), the vehicle state equation can be described as(4)
v
˙
y
γ
˙
=
-
K
f
+
K
r
m
v
x
-
a
K
f
-
b
K
r
m
v
x
-
v
x
-
a
K
f
-
b
K
r
I
z
v
x
-
a
2
K
f
+
b
2
K
r
I
z
v
x
v
y
γ
+
K
f
m
a
K
f
I
z
δ
.In practice, only the lateral accelerationa
y and yaw rate γ are available by vehicle sensors, but the measurable outputs cannot be expressed by the linear combination of state variables v
y and γ. In fact, the measurable output a
y is equal to the sum of vehicle speed v
x multiplied by yaw rate γ and the derivative of lateral velocity v
˙
y (a
y
=
v
x
γ
+
v
˙
y). Therefore, it is necessary to extend (4) so as to obtain the state variable v
˙
y. Taking the derivative of (4), an extended state equation, including the derivative of lateral velocity, can be obtained as follows: (5)
v
¨
y
v
˙
y
γ
¨
γ
˙
=
-
K
f
+
K
r
m
v
x
0
-
a
K
f
-
b
K
r
m
v
x
-
v
x
0
0
-
K
f
+
K
r
m
v
x
0
-
a
K
f
-
b
K
r
m
v
x
-
v
x
-
a
K
f
-
b
K
r
I
z
U
0
-
a
2
K
f
+
b
2
K
r
I
z
v
x
0
0
-
a
K
f
-
b
K
r
I
z
v
x
0
-
a
2
K
f
+
b
2
K
r
I
z
v
x
v
˙
y
v
y
γ
˙
γ
+
K
f
m
0
0
K
f
m
a
K
f
I
z
0
0
a
K
f
I
z
δ
˙
δ
.Given a longitudinal speedv
x, the vehicle lateral dynamics model is a linear time-invariant system. When the measured steering angle has certain error or noise, the derivative may be very noisy. Therefore, the steering angle signal should be filtered by a mid-value filtering smoother before calculating the steering angular velocity. However, assumption of LTI vehicle model cannot be strictly ensured in practical application. Vehicle speed is an average value, which fluctuates slightly over time in the road tests. Moreover, the cornering stiffness is time-varying in strict sense. As seen in Figure 2, the tire cornering stiffness changes slowly with the transfer of vertical load or the increase of tire slip angle even in the linear domain.Figure 2
The relationship between tire lateral force and slip angle (F
z is the vertical load).In consideration of these time-varying parameters, a linear time-varying vehicle model is more reasonable and accurate to describe the vehicle performance. As a result, the vehicle model can be expressed as the following form:(6)
x
˙
=
G
t
x
+
H
t
u
,
y
=
C
t
x
+
D
t
u
,where x
=
[
v
˙
y
v
y
γ
˙
γ
]
T, u
=
[
δ
˙
δ
]
T, y
=
[
a
y
γ
]
T, and system matrices are(7)
G
t
=
-
K
f
t
+
K
r
t
m
v
x
t
0
-
a
K
f
t
-
b
K
r
t
m
v
x
t
-
v
x
t
0
0
-
K
f
t
+
K
r
t
m
v
x
t
0
-
a
K
f
t
-
b
K
r
t
m
v
x
t
-
v
x
t
-
a
K
f
t
-
b
K
r
t
I
z
v
x
t
0
-
a
2
K
f
t
+
b
2
K
r
t
I
z
v
x
t
0
0
-
a
K
f
t
-
b
K
r
t
I
z
v
x
t
0
-
a
2
K
f
t
+
b
2
K
r
t
I
z
v
x
t
,
H
t
=
K
f
t
m
0
a
K
f
t
I
z
0
0
K
f
t
m
0
a
K
f
t
I
z
T
,
C
t
=
1
0
0
v
x
t
0
0
0
1
,
D
t
=
0
.Considering the time-varying characteristics of vehicle speed and cornering stiffness, this model structure can describe vehicle behavior accurately and the application scope is not strictly limited to the linear range for lateral acceleration.One precondition of theR
P
B
S
I
D
o
p
t method is that the system to be identified should be minimal, and this means system (6) should be observable and controllable for any time t. However, the accurate system matrices are not known priorly, and they need to be estimated and updated by input and output data. Therefore, the proof of minimal system is hard to be given in strict sense. In fact, the time-varying parameters K
f
(
t
), K
r
(
t
), and v
x
(
t
) in system (6) change slowly near their nominal values in standard road tests, and these parameters are also bounded. Thus, the precondition is approximately demonstrated when parameters are at nominal values.Define the observability matrixΓ and controllability matrix Δ as(8)
Γ
=
C
T
G
T
C
T
G
2
T
C
T
G
3
T
C
T
T
,
(9)
Δ
=
H
G
H
G
2
H
G
3
H
.The specific proof of minimal system is given in the Appendix. The system is observable and controllable, and this conclusion is confirmed by the fact that steering wheel angle and angular velocity are the major inputs to keep directional control steady. To handle the black-box problem of the identifiable 4th-order vehicle model, the previous identification methods are not suitable. In next subsection, a recursive subspace identification method is introduced.


### 2.2. Review ofR
P
B
S
I
D
o
p
t Algorithm

TheP
B
S
I
D
o
p
t method introduced in [13, 22] is becoming an effective and well-established technique for real LTI systems [23, 24]. In order to avoid the heavy computation of SVD, OPAST is employed to update the extended observability matrix and implement the recursive P
B
S
I
D
o
p
t (R
P
B
S
I
D
o
p
t) identification.Considering the process and measurement noises, the discrete form of an LTI system given can be expressed as(10)
x
t
+
1
=
A
x
t
+
B
u
t
+
K
e
t
,
(11)
y
t
=
C
x
t
+
D
u
t
+
e
t
,where u
(
t
)
∈
R
m and y
(
t
)
∈
R
l are input and output sequences, respectively. x
(
t
)
∈
R
n is a state vector, and e
(
t
)
∈
R
l is defined as e
(
t
)
=
y
(
t
)
-
y
^
(
t
∣
t
-
1
), where y
^
(
t
∣
t
-
1
) means the estimation of y
(
t
) by using the data before moment t
-
1. For the proposed vehicle structure, the matrix D is equal to 0. In addition, the following assumptions for the system in (10) and (11) hold.(A1)
The system is minimal, that is, the dimensionn of x
(
t
) is as small as possible.
(A2)
The matrixA
-
K
C is asymptotically stable.
(A3)
The input signalu
(
t
) is persistently exciting (PE).The specific definition of PE condition can be referred to in literature [25]. The judgment of PE condition mainly depends on the input spectrum, and the input signal should contain sufficiently distinct frequencies. Furthermore, a discriminated method for the finite-dimension input is proposed by Katayama [26]. In fact, the PE condition of vehicle steering input has been researched in our previous study [9], and it is satisfied for most standard road tests. In fact, the conditions should guarantee (a) a proper recovery of the Markov parameters from (26) and (b) strict convexity of the objective function in (34). In a batch identification, the matrix z
h should have full row-rank to achieve (a). And in recursive mode, the condition (a) for convergence of the RLS is still related to the sequence z
h but is a bit stronger [27]. However, these considerations are too theoretical, and the PE condition is not the research emphasis in this paper. So the above conditions (a) and (b) are assumed to be satisfied, and all signals have the necessary properties allowing for the proper recovery of the different matrices.Accordingly, the key idea of the vehicle model identification is then transferred to determination of the system matricesA,B,C,K, and the innovation process e
(
t
) by the measurements of input and output (I/O) sequences u
(
k
)
,
y
(
k
)
k
=
1
k
=
t.The system expressed by (10) and (11) can be rewritten in a predictor form as(12)
x
t
+
1
=
A
-
x
t
+
B
-
z
t
,
(13)
y
t
=
C
x
t
+
e
t
,where z
t
:
=
[
u
T
(
t
)
y
T
(
t
)
]
T
∈
R
4 and A
-
:
=
A
-
K
C, B
-
:
=
[
B
K
].AnR
P
B
S
I
D
o
p
t algorithm is employed to handle the vehicle model identification problem. The specific derivation process of the R
P
B
S
I
D
o
p
t algorithm can be referred to in the literature [13, 20]. The 3-step calculation process is given below.Repeatedly substituting the state equation in (12), it can be derived as(14)
x
t
=
A
-
h
x
t
-
h
+
∑
i
=
1
h
A
-
i
-
1
B
-
z
t
-
i
.According to the assumption (A2), A
-
h approaches to0 when the past horizon h is sufficiently large. This means that the first term on right side of (14) is negligible. Hence, (14) can be expressed as the following form:(15)
x
t
≈
∑
i
=
1
h
A
-
i
-
1
B
-
z
t
-
i
≔
L
h
z
h
t
,where L
h
=
[
A
-
h
-
1
B
-
A
-
h
-
2
B
-
⋯
B
-
] and z
h
t
=
[
z
T
(
t
-
h
)
z
T
(
t
-
h
+
1
)
⋯
z
T
(
t
-
1
)
]
T.Iteration of (12) and (13) results in the subspace matrix formulation of the following form as(16)
Y
f
=
Γ
-
f
X
t
+
H
-
f
Z
f
-
1
+
E
f
,where X
t
=
[
x
(
t
+
1
)
x
(
t
+
2
)
⋯
x
(
t
+
i
)
], Y
f, Z
f
-
1, and E
f are block-Hankel matrices, and(17)
Y
f
=
y
t
+
1
y
t
+
2
⋯
y
t
+
i
y
t
+
2
y
t
+
3
⋯
y
t
+
i
+
1
⋮
⋮
⋱
⋮
y
t
+
f
y
t
+
f
+
1
⋯
y
t
+
f
+
i
-
1
,
Z
f
-
1
=
z
t
+
1
z
t
+
2
⋯
z
t
+
i
z
t
+
2
z
t
+
3
⋯
z
t
+
i
+
1
⋮
⋮
⋱
⋮
z
t
+
f
-
1
z
t
+
f
⋯
z
t
+
f
+
i
-
2
,
E
f
=
e
t
+
1
e
t
+
2
⋯
e
t
+
i
e
t
+
2
e
t
+
3
⋯
e
t
+
i
+
1
⋮
⋮
⋱
⋮
e
t
+
f
e
t
+
f
+
1
⋯
e
t
+
f
+
i
-
1
.The extended observability matrix is given as Γ
-
f
=
[
C
T
(
C
A
-
)
T
⋯
(
C
A
-
f
-
1
)
T
]
T, and the lower triangular block-Toeplitz matrix H
-
f is given by(18)
H
-
f
=
0
0
⋯
0
C
B
-
0
⋯
0
C
A
-
B
-
C
B
-
⋱
0
⋮
⋮
⋱
⋮
C
A
-
f
-
2
B
-
C
A
-
f
-
3
B
-
⋯
C
B
-
.According to (15), the state matrix X
t can be expressed as(19)
X
t
=
L
h
Z
h
,where Z
h
=
[
z
h
(
t
)
z
h
(
t
+
1
)
⋯
z
h
(
t
+
i
-
1
)
].Substituting (19) into (16), then(20)
Y
f
=
Γ
-
f
L
h
Z
h
+
H
-
f
Z
f
-
1
+
E
f
.The matrix M is defined by(21)
M
≔
Γ
-
f
L
h
=
C
C
A
-
⋮
C
A
-
f
-
1
A
-
h
-
1
B
-
A
-
h
-
2
B
-
⋯
B
-
=
m
h
m
h
-
1
⋯
m
1
m
h
+
1
m
h
⋯
m
2
⋮
⋮
⋱
⋮
m
h
+
f
-
1
m
h
+
f
-
2
⋯
m
f
,where m
i
:
=
C
A
-
i
-
1
B
-  (
i
=
1,2
,
…
) is called the predictor Markov parameters. As A
-
h
=
0 whenever h is sufficiently large, m
i
=
0  (
i
>
h
). Then, the matrix M in (21) can be simplified to the following form:(22)
M
=
m
h
m
h
-
1
⋯
m
1
0
m
h
⋯
m
2
⋮
⋮
⋱
⋮
0
0
⋯
m
f
.According to (16) and (20), an important identity for P
B
S
I
D
o
p
t algorithm can be found as below:(23)
Γ
-
f
X
t
=
M
Z
h
.And the first column vector form of (23) is(24)
Γ
-
f
x
t
+
1
=
M
z
h
t
.The identity is the key equation for P
B
S
I
D
o
p
t algorithm.Step 1 (the estimation of Markov parametersm
i).
Substitution of (15) into (13) results in the following form:(25)
y
t
=
∑
i
=
1
h
C
A
-
i
-
1
B
-
z
t
-
i
+
e
t
=
∑
i
=
1
h
m
i
z
t
-
i
+
e
t
.According to (13), the output vector y
(
t
) can be expressed as(26)
y
t
=
Θ
z
h
t
+
e
t
,where Θ
≔
m
h
m
h
-
1
⋯
m
1, and the element m
i
:
=
C
A
-
i
-
1
B
-  i
=
1
,
…
,
h is called predictor Markov parameters: z
h
t
=
[
z
T
(
t
-
h
)
z
T
(
t
-
h
+
1
)
⋯
z
T
(
t
-
1
)
]
T.
The Markov parametersm
i can be estimated by RLS (Recursive Least Square) filter with forgetting factor λ
1.Step 2 (the update of the extended observability matrixΓ
-
f).
An OPAST method is proposed to update the extended matrixΓ
-
f. As discussed in [21], it can guarantee the orthonormality of the weighted matrix and converges globally.
Consider the unconstrained criteria with a random vectorξ
∈
R
m as follows:(27)
J
W
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
-
W
t
W
T
t
ξ
k
2
,where λ
2 is a forgetting factor and W
∈
R
m
×
n  (
m
>
n
) is a full rank matrix.Theorem 1 (see [19]).
The global minimumJ
(
W
) is attained if and only if W
=
Q
-
T, where Q
- contains the n dominating eigenvectors of R
ξ. Here T is an arbitrary unitary matrix. Furthermore, all other stationary points are saddle points.The above theorem is also applicable ifR
ξ is replaced with(28)
R
^
ξ
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
ξ
T
k
.The key idea of PAST is to replaceW
T
(
t
)
ξ
(
k
) in (27) with(29)
h
k
=
W
T
k
-
1
x
k
.The projection approximation is valid if the weighted matrix is slowly varying. Then the criterion of quadratic form with respect to W
(
t
) is rewritten as follows:(30)
J
-
W
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
-
W
t
h
k
2
,which is minimized by(31)
W
t
=
R
^
ξ
h
t
R
^
h
-
1
t
,where R
^
ξ
h
(
t
)
=
∑
k
=
1
t
λ
2
t
-
k
ξ
(
k
)
h
T
(
k
) and R
^
h
(
t
)
=
∑
k
=
1
t
λ
2
t
-
k
h
(
k
)
h
T
(
k
) are covariance matrices, and it is assumed that the inverse of matrix R
^
h
(
t
) exists.Since the criterionJ
-
(
W
) is used for approximation of J
(
W
), the estimation may be oscillating without convergence in some cases. For alleviation of the divergence, the OPAST method adds an orthonormalization step at each iteration:(32)
W
t
=
W
t
W
T
t
W
t
-
1
/
2
.As observed in identity equation (24), the column space of M
z
h
(
t
) is consistent with the column space of Γ
-
f. M
z
h
(
t
) is fed into objective function (27) as the vector ξ, and then the output W
(
k
) is corresponding to the extended observability matrix Γ
-
f. And this rationale is widely used in previous literature [17, 28] for updating observability matrix Γ
-
f.Step 3 (the estimation of the system matricesA
,
B
,
C
,
K).
Calculate the estimated matrixC
^. It can be extracted directly and equals the front l rows of Γ
-
f according to the definition:(33)
C
^
=
Γ
-
f
1
:
l
,
:
,where Γ
-
f
(
1
:
l
,
:
) is a MATLAB expression, and it indicates the front l rows.
According to (10), the system matrices {
A
^
(
t
)
B
^
(
t
)
K
^
(
t
)
} can be estimated by minimizing the following criterion:(34)
min
⁡
J
=
∑
k
=
1
t
x
k
-
A
t
B
t
K
t
x
k
-
1
u
k
-
1
e
k
-
1
2
,where the innovation e
(
k
)
=
y
(
k
)
-
C
^
(
k
)
x
(
k
). And the states x can be estimated according to (24), where matrices Γ
-
f and M are obtained before and the operation “†” donates Moore-Penrose inverse.


### 2.3. Vehicle Model Identification Process

Due to the time-varying parameters in (6) changing slowly near their nominal values under standard road tests, the vehicle model is approximately regarded as LTI system at each moment. And the coefficient matrices in discrete-time system are A
=
e
G
T and B
=
∫
0
T
e
G
t
H
d
t, where T is a sampling interval. Through above approximation, the R
P
B
S
I
D
o
p
t method can be applied to estimate and update the model matrices at each moment t, and the detailed process is given in Figure 3. As shown, the currently available input and output of real vehicle system are collected and imported into the identification module. By the 3-step R
P
B
S
I
D
o
p
t calculation process, the vehicle model matrices, used to predict the vehicle outputs, can be identified in real-time.Figure 3
Vehicle model identification process.


### 2.4. Vehicle Road Tests

In order to verify the practicability and effectiveness of the identification model, four vehicle road tests are carried out, including double lane change test (ISO 3888-1:1999), pulse input test (ISO7401), pylon course slalom test (Chinese standard GB/T 6323.1-94), and step input test (ISO7401). As shown in Figure4, steering wheel angle δ is measured by a gauge as one of the input signals, and the steering angular velocity δ
˙ needs to be calculated indirectly by δ and sampling time. It is noted that before calculating the steering angular velocity, the steering angle signal should be filtered by smoother in order to relieve the divergence caused by certain error or noise. Yaw rate γ is measured by gyroscope as one of the output signals. The lateral acceleration a
y can be obtained by accelerometer as another output signal, which is the sum of vehicle speed multiplied by yaw rate and the derivative of lateral velocity. In order to obtain accurate measurements, the accelerometer and gyroscope are placed on the vehicle center of gravity. It should be noted that the vertical gyro “VG400CD-200” by Moog Crossbow is a sensor cluster of accelerometer and gyroscope. And a low-pass Kaiser Filter is used for the collected signals. During the road tests, sample frequency is set as 100 Hz. And Table 1 lists the characteristic of the instruments in the road tests.Table 1
List of the characteristics of instruments in the road test.
Sensor
Type
Range
Accuracy
GPS
Dewe-VGPS-200
—
0.1 m/s
Data acquisition instrument
Dewe-501
—
—
Steering wheel gauge
MSW-1
±720 deg
0.1 deg
Vertical gyroscope
Moog-VG400CD-200
Acceleration: ±4 g
1 mg
Angular rate: ±200 deg/s
0.02 deg/sFigure 4
The vehicle and instruments in road test.


## 2.1. Identifiable Vehicle Model Structure

A commonly used model in vehicle lateral dynamics research is the single track model, which can simulate the lateral behavior under many conditions. Under the assumptions of small steering angle and same slip angles of left and right tires, and ignoring the roll motion, the differential equations of vehicle lateral motion and yaw motion are expressed as below:(1)
m
v
˙
y
+
v
x
γ
=
F
f
+
F
r
,
I
z
γ
˙
=
a
F
f
-
b
F
r
,where m is the vehicle mass, v
x is the longitudinal speed, v
y is the lateral speed of the CG of the vehicle, γ is the yaw rate, I
z is the vehicle yaw moment inertia, and F
f and F
r are front and rear tire lateral forces (Figure 1).Figure 1
The 2-DOF single track model.For linear tire model under small tire slip angle, the lateral forces can be written as(2)
F
f
=
-
K
f
β
f
,
F
r
=
-
K
r
β
r
,where K
f and K
r are the average cornering stiffness of front and rear axles and β
f and β
r are the front and rear tire slip angles, respectively. The tire slip angles can be expressed approximatively as(3)
β
f
=
-
δ
+
v
y
+
a
γ
v
x
,
β
r
=
v
y
-
b
γ
v
x
.Substituting (2) and (3) into (1), the vehicle state equation can be described as(4)
v
˙
y
γ
˙
=
-
K
f
+
K
r
m
v
x
-
a
K
f
-
b
K
r
m
v
x
-
v
x
-
a
K
f
-
b
K
r
I
z
v
x
-
a
2
K
f
+
b
2
K
r
I
z
v
x
v
y
γ
+
K
f
m
a
K
f
I
z
δ
.In practice, only the lateral accelerationa
y and yaw rate γ are available by vehicle sensors, but the measurable outputs cannot be expressed by the linear combination of state variables v
y and γ. In fact, the measurable output a
y is equal to the sum of vehicle speed v
x multiplied by yaw rate γ and the derivative of lateral velocity v
˙
y (a
y
=
v
x
γ
+
v
˙
y). Therefore, it is necessary to extend (4) so as to obtain the state variable v
˙
y. Taking the derivative of (4), an extended state equation, including the derivative of lateral velocity, can be obtained as follows: (5)
v
¨
y
v
˙
y
γ
¨
γ
˙
=
-
K
f
+
K
r
m
v
x
0
-
a
K
f
-
b
K
r
m
v
x
-
v
x
0
0
-
K
f
+
K
r
m
v
x
0
-
a
K
f
-
b
K
r
m
v
x
-
v
x
-
a
K
f
-
b
K
r
I
z
U
0
-
a
2
K
f
+
b
2
K
r
I
z
v
x
0
0
-
a
K
f
-
b
K
r
I
z
v
x
0
-
a
2
K
f
+
b
2
K
r
I
z
v
x
v
˙
y
v
y
γ
˙
γ
+
K
f
m
0
0
K
f
m
a
K
f
I
z
0
0
a
K
f
I
z
δ
˙
δ
.Given a longitudinal speedv
x, the vehicle lateral dynamics model is a linear time-invariant system. When the measured steering angle has certain error or noise, the derivative may be very noisy. Therefore, the steering angle signal should be filtered by a mid-value filtering smoother before calculating the steering angular velocity. However, assumption of LTI vehicle model cannot be strictly ensured in practical application. Vehicle speed is an average value, which fluctuates slightly over time in the road tests. Moreover, the cornering stiffness is time-varying in strict sense. As seen in Figure 2, the tire cornering stiffness changes slowly with the transfer of vertical load or the increase of tire slip angle even in the linear domain.Figure 2
The relationship between tire lateral force and slip angle (F
z is the vertical load).In consideration of these time-varying parameters, a linear time-varying vehicle model is more reasonable and accurate to describe the vehicle performance. As a result, the vehicle model can be expressed as the following form:(6)
x
˙
=
G
t
x
+
H
t
u
,
y
=
C
t
x
+
D
t
u
,where x
=
[
v
˙
y
v
y
γ
˙
γ
]
T, u
=
[
δ
˙
δ
]
T, y
=
[
a
y
γ
]
T, and system matrices are(7)
G
t
=
-
K
f
t
+
K
r
t
m
v
x
t
0
-
a
K
f
t
-
b
K
r
t
m
v
x
t
-
v
x
t
0
0
-
K
f
t
+
K
r
t
m
v
x
t
0
-
a
K
f
t
-
b
K
r
t
m
v
x
t
-
v
x
t
-
a
K
f
t
-
b
K
r
t
I
z
v
x
t
0
-
a
2
K
f
t
+
b
2
K
r
t
I
z
v
x
t
0
0
-
a
K
f
t
-
b
K
r
t
I
z
v
x
t
0
-
a
2
K
f
t
+
b
2
K
r
t
I
z
v
x
t
,
H
t
=
K
f
t
m
0
a
K
f
t
I
z
0
0
K
f
t
m
0
a
K
f
t
I
z
T
,
C
t
=
1
0
0
v
x
t
0
0
0
1
,
D
t
=
0
.Considering the time-varying characteristics of vehicle speed and cornering stiffness, this model structure can describe vehicle behavior accurately and the application scope is not strictly limited to the linear range for lateral acceleration.One precondition of theR
P
B
S
I
D
o
p
t method is that the system to be identified should be minimal, and this means system (6) should be observable and controllable for any time t. However, the accurate system matrices are not known priorly, and they need to be estimated and updated by input and output data. Therefore, the proof of minimal system is hard to be given in strict sense. In fact, the time-varying parameters K
f
(
t
), K
r
(
t
), and v
x
(
t
) in system (6) change slowly near their nominal values in standard road tests, and these parameters are also bounded. Thus, the precondition is approximately demonstrated when parameters are at nominal values.Define the observability matrixΓ and controllability matrix Δ as(8)
Γ
=
C
T
G
T
C
T
G
2
T
C
T
G
3
T
C
T
T
,
(9)
Δ
=
H
G
H
G
2
H
G
3
H
.The specific proof of minimal system is given in the Appendix. The system is observable and controllable, and this conclusion is confirmed by the fact that steering wheel angle and angular velocity are the major inputs to keep directional control steady. To handle the black-box problem of the identifiable 4th-order vehicle model, the previous identification methods are not suitable. In next subsection, a recursive subspace identification method is introduced.


## 2.2. Review ofR
P
B
S
I
D
o
p
t Algorithm

TheP
B
S
I
D
o
p
t method introduced in [13, 22] is becoming an effective and well-established technique for real LTI systems [23, 24]. In order to avoid the heavy computation of SVD, OPAST is employed to update the extended observability matrix and implement the recursive P
B
S
I
D
o
p
t (R
P
B
S
I
D
o
p
t) identification.Considering the process and measurement noises, the discrete form of an LTI system given can be expressed as(10)
x
t
+
1
=
A
x
t
+
B
u
t
+
K
e
t
,
(11)
y
t
=
C
x
t
+
D
u
t
+
e
t
,where u
(
t
)
∈
R
m and y
(
t
)
∈
R
l are input and output sequences, respectively. x
(
t
)
∈
R
n is a state vector, and e
(
t
)
∈
R
l is defined as e
(
t
)
=
y
(
t
)
-
y
^
(
t
∣
t
-
1
), where y
^
(
t
∣
t
-
1
) means the estimation of y
(
t
) by using the data before moment t
-
1. For the proposed vehicle structure, the matrix D is equal to 0. In addition, the following assumptions for the system in (10) and (11) hold.(A1)
The system is minimal, that is, the dimensionn of x
(
t
) is as small as possible.
(A2)
The matrixA
-
K
C is asymptotically stable.
(A3)
The input signalu
(
t
) is persistently exciting (PE).The specific definition of PE condition can be referred to in literature [25]. The judgment of PE condition mainly depends on the input spectrum, and the input signal should contain sufficiently distinct frequencies. Furthermore, a discriminated method for the finite-dimension input is proposed by Katayama [26]. In fact, the PE condition of vehicle steering input has been researched in our previous study [9], and it is satisfied for most standard road tests. In fact, the conditions should guarantee (a) a proper recovery of the Markov parameters from (26) and (b) strict convexity of the objective function in (34). In a batch identification, the matrix z
h should have full row-rank to achieve (a). And in recursive mode, the condition (a) for convergence of the RLS is still related to the sequence z
h but is a bit stronger [27]. However, these considerations are too theoretical, and the PE condition is not the research emphasis in this paper. So the above conditions (a) and (b) are assumed to be satisfied, and all signals have the necessary properties allowing for the proper recovery of the different matrices.Accordingly, the key idea of the vehicle model identification is then transferred to determination of the system matricesA,B,C,K, and the innovation process e
(
t
) by the measurements of input and output (I/O) sequences u
(
k
)
,
y
(
k
)
k
=
1
k
=
t.The system expressed by (10) and (11) can be rewritten in a predictor form as(12)
x
t
+
1
=
A
-
x
t
+
B
-
z
t
,
(13)
y
t
=
C
x
t
+
e
t
,where z
t
:
=
[
u
T
(
t
)
y
T
(
t
)
]
T
∈
R
4 and A
-
:
=
A
-
K
C, B
-
:
=
[
B
K
].AnR
P
B
S
I
D
o
p
t algorithm is employed to handle the vehicle model identification problem. The specific derivation process of the R
P
B
S
I
D
o
p
t algorithm can be referred to in the literature [13, 20]. The 3-step calculation process is given below.Repeatedly substituting the state equation in (12), it can be derived as(14)
x
t
=
A
-
h
x
t
-
h
+
∑
i
=
1
h
A
-
i
-
1
B
-
z
t
-
i
.According to the assumption (A2), A
-
h approaches to0 when the past horizon h is sufficiently large. This means that the first term on right side of (14) is negligible. Hence, (14) can be expressed as the following form:(15)
x
t
≈
∑
i
=
1
h
A
-
i
-
1
B
-
z
t
-
i
≔
L
h
z
h
t
,where L
h
=
[
A
-
h
-
1
B
-
A
-
h
-
2
B
-
⋯
B
-
] and z
h
t
=
[
z
T
(
t
-
h
)
z
T
(
t
-
h
+
1
)
⋯
z
T
(
t
-
1
)
]
T.Iteration of (12) and (13) results in the subspace matrix formulation of the following form as(16)
Y
f
=
Γ
-
f
X
t
+
H
-
f
Z
f
-
1
+
E
f
,where X
t
=
[
x
(
t
+
1
)
x
(
t
+
2
)
⋯
x
(
t
+
i
)
], Y
f, Z
f
-
1, and E
f are block-Hankel matrices, and(17)
Y
f
=
y
t
+
1
y
t
+
2
⋯
y
t
+
i
y
t
+
2
y
t
+
3
⋯
y
t
+
i
+
1
⋮
⋮
⋱
⋮
y
t
+
f
y
t
+
f
+
1
⋯
y
t
+
f
+
i
-
1
,
Z
f
-
1
=
z
t
+
1
z
t
+
2
⋯
z
t
+
i
z
t
+
2
z
t
+
3
⋯
z
t
+
i
+
1
⋮
⋮
⋱
⋮
z
t
+
f
-
1
z
t
+
f
⋯
z
t
+
f
+
i
-
2
,
E
f
=
e
t
+
1
e
t
+
2
⋯
e
t
+
i
e
t
+
2
e
t
+
3
⋯
e
t
+
i
+
1
⋮
⋮
⋱
⋮
e
t
+
f
e
t
+
f
+
1
⋯
e
t
+
f
+
i
-
1
.The extended observability matrix is given as Γ
-
f
=
[
C
T
(
C
A
-
)
T
⋯
(
C
A
-
f
-
1
)
T
]
T, and the lower triangular block-Toeplitz matrix H
-
f is given by(18)
H
-
f
=
0
0
⋯
0
C
B
-
0
⋯
0
C
A
-
B
-
C
B
-
⋱
0
⋮
⋮
⋱
⋮
C
A
-
f
-
2
B
-
C
A
-
f
-
3
B
-
⋯
C
B
-
.According to (15), the state matrix X
t can be expressed as(19)
X
t
=
L
h
Z
h
,where Z
h
=
[
z
h
(
t
)
z
h
(
t
+
1
)
⋯
z
h
(
t
+
i
-
1
)
].Substituting (19) into (16), then(20)
Y
f
=
Γ
-
f
L
h
Z
h
+
H
-
f
Z
f
-
1
+
E
f
.The matrix M is defined by(21)
M
≔
Γ
-
f
L
h
=
C
C
A
-
⋮
C
A
-
f
-
1
A
-
h
-
1
B
-
A
-
h
-
2
B
-
⋯
B
-
=
m
h
m
h
-
1
⋯
m
1
m
h
+
1
m
h
⋯
m
2
⋮
⋮
⋱
⋮
m
h
+
f
-
1
m
h
+
f
-
2
⋯
m
f
,where m
i
:
=
C
A
-
i
-
1
B
-  (
i
=
1,2
,
…
) is called the predictor Markov parameters. As A
-
h
=
0 whenever h is sufficiently large, m
i
=
0  (
i
>
h
). Then, the matrix M in (21) can be simplified to the following form:(22)
M
=
m
h
m
h
-
1
⋯
m
1
0
m
h
⋯
m
2
⋮
⋮
⋱
⋮
0
0
⋯
m
f
.According to (16) and (20), an important identity for P
B
S
I
D
o
p
t algorithm can be found as below:(23)
Γ
-
f
X
t
=
M
Z
h
.And the first column vector form of (23) is(24)
Γ
-
f
x
t
+
1
=
M
z
h
t
.The identity is the key equation for P
B
S
I
D
o
p
t algorithm.Step 1 (the estimation of Markov parametersm
i).
Substitution of (15) into (13) results in the following form:(25)
y
t
=
∑
i
=
1
h
C
A
-
i
-
1
B
-
z
t
-
i
+
e
t
=
∑
i
=
1
h
m
i
z
t
-
i
+
e
t
.According to (13), the output vector y
(
t
) can be expressed as(26)
y
t
=
Θ
z
h
t
+
e
t
,where Θ
≔
m
h
m
h
-
1
⋯
m
1, and the element m
i
:
=
C
A
-
i
-
1
B
-  i
=
1
,
…
,
h is called predictor Markov parameters: z
h
t
=
[
z
T
(
t
-
h
)
z
T
(
t
-
h
+
1
)
⋯
z
T
(
t
-
1
)
]
T.
The Markov parametersm
i can be estimated by RLS (Recursive Least Square) filter with forgetting factor λ
1.Step 2 (the update of the extended observability matrixΓ
-
f).
An OPAST method is proposed to update the extended matrixΓ
-
f. As discussed in [21], it can guarantee the orthonormality of the weighted matrix and converges globally.
Consider the unconstrained criteria with a random vectorξ
∈
R
m as follows:(27)
J
W
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
-
W
t
W
T
t
ξ
k
2
,where λ
2 is a forgetting factor and W
∈
R
m
×
n  (
m
>
n
) is a full rank matrix.Theorem 1 (see [19]).
The global minimumJ
(
W
) is attained if and only if W
=
Q
-
T, where Q
- contains the n dominating eigenvectors of R
ξ. Here T is an arbitrary unitary matrix. Furthermore, all other stationary points are saddle points.The above theorem is also applicable ifR
ξ is replaced with(28)
R
^
ξ
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
ξ
T
k
.The key idea of PAST is to replaceW
T
(
t
)
ξ
(
k
) in (27) with(29)
h
k
=
W
T
k
-
1
x
k
.The projection approximation is valid if the weighted matrix is slowly varying. Then the criterion of quadratic form with respect to W
(
t
) is rewritten as follows:(30)
J
-
W
t
=
∑
k
=
1
t
λ
2
t
-
k
ξ
k
-
W
t
h
k
2
,which is minimized by(31)
W
t
=
R
^
ξ
h
t
R
^
h
-
1
t
,where R
^
ξ
h
(
t
)
=
∑
k
=
1
t
λ
2
t
-
k
ξ
(
k
)
h
T
(
k
) and R
^
h
(
t
)
=
∑
k
=
1
t
λ
2
t
-
k
h
(
k
)
h
T
(
k
) are covariance matrices, and it is assumed that the inverse of matrix R
^
h
(
t
) exists.Since the criterionJ
-
(
W
) is used for approximation of J
(
W
), the estimation may be oscillating without convergence in some cases. For alleviation of the divergence, the OPAST method adds an orthonormalization step at each iteration:(32)
W
t
=
W
t
W
T
t
W
t
-
1
/
2
.As observed in identity equation (24), the column space of M
z
h
(
t
) is consistent with the column space of Γ
-
f. M
z
h
(
t
) is fed into objective function (27) as the vector ξ, and then the output W
(
k
) is corresponding to the extended observability matrix Γ
-
f. And this rationale is widely used in previous literature [17, 28] for updating observability matrix Γ
-
f.Step 3 (the estimation of the system matricesA
,
B
,
C
,
K).
Calculate the estimated matrixC
^. It can be extracted directly and equals the front l rows of Γ
-
f according to the definition:(33)
C
^
=
Γ
-
f
1
:
l
,
:
,where Γ
-
f
(
1
:
l
,
:
) is a MATLAB expression, and it indicates the front l rows.
According to (10), the system matrices {
A
^
(
t
)
B
^
(
t
)
K
^
(
t
)
} can be estimated by minimizing the following criterion:(34)
min
⁡
J
=
∑
k
=
1
t
x
k
-
A
t
B
t
K
t
x
k
-
1
u
k
-
1
e
k
-
1
2
,where the innovation e
(
k
)
=
y
(
k
)
-
C
^
(
k
)
x
(
k
). And the states x can be estimated according to (24), where matrices Γ
-
f and M are obtained before and the operation “†” donates Moore-Penrose inverse.


## 2.3. Vehicle Model Identification Process

Due to the time-varying parameters in (6) changing slowly near their nominal values under standard road tests, the vehicle model is approximately regarded as LTI system at each moment. And the coefficient matrices in discrete-time system are A
=
e
G
T and B
=
∫
0
T
e
G
t
H
d
t, where T is a sampling interval. Through above approximation, the R
P
B
S
I
D
o
p
t method can be applied to estimate and update the model matrices at each moment t, and the detailed process is given in Figure 3. As shown, the currently available input and output of real vehicle system are collected and imported into the identification module. By the 3-step R
P
B
S
I
D
o
p
t calculation process, the vehicle model matrices, used to predict the vehicle outputs, can be identified in real-time.Figure 3
Vehicle model identification process.


## 2.4. Vehicle Road Tests

In order to verify the practicability and effectiveness of the identification model, four vehicle road tests are carried out, including double lane change test (ISO 3888-1:1999), pulse input test (ISO7401), pylon course slalom test (Chinese standard GB/T 6323.1-94), and step input test (ISO7401). As shown in Figure4, steering wheel angle δ is measured by a gauge as one of the input signals, and the steering angular velocity δ
˙ needs to be calculated indirectly by δ and sampling time. It is noted that before calculating the steering angular velocity, the steering angle signal should be filtered by smoother in order to relieve the divergence caused by certain error or noise. Yaw rate γ is measured by gyroscope as one of the output signals. The lateral acceleration a
y can be obtained by accelerometer as another output signal, which is the sum of vehicle speed multiplied by yaw rate and the derivative of lateral velocity. In order to obtain accurate measurements, the accelerometer and gyroscope are placed on the vehicle center of gravity. It should be noted that the vertical gyro “VG400CD-200” by Moog Crossbow is a sensor cluster of accelerometer and gyroscope. And a low-pass Kaiser Filter is used for the collected signals. During the road tests, sample frequency is set as 100 Hz. And Table 1 lists the characteristic of the instruments in the road tests.Table 1
List of the characteristics of instruments in the road test.
Sensor
Type
Range
Accuracy
GPS
Dewe-VGPS-200
—
0.1 m/s
Data acquisition instrument
Dewe-501
—
—
Steering wheel gauge
MSW-1
±720 deg
0.1 deg
Vertical gyroscope
Moog-VG400CD-200
Acceleration: ±4 g
1 mg
Angular rate: ±200 deg/s
0.02 deg/sFigure 4
The vehicle and instruments in road test.


## 3. Vehicle Model Identification and Test Validation

### 3.1. The Choice of Algorithm Parameters

Before vehicle identification, there are four constant parameters in the algorithm to be determined, and they are the past and future horizonsh and f and forgetting factors λ
1 and λ
2, respectively. These parameters play important roles in the convergence rate and identification precision. Based on vehicle pulse input test, the influences of the pending parameters are studied, and it could guide for an appropriate choice of these parameters.Generally, the past and future horizons are user choices in the subspace algorithm. And the future horizonf is less than or equal to the past horizon h (f
≤
h). As pointed out by Houtzager et al. [29], the truncation error decreased with the increase of past horizon h, but the computational time became larger. On the other hand, Chiuso [30] suggested that a large future horizon f was often better for identification when the signal-to-noise ratio was low. And Houtzager et al. [29] also indicated that when premultiplied weighting matrix was not applied for R
P
B
S
I
D
o
p
t algorithm, the high identification accuracy was obtained when the future horizon was close to the past horizon. Therefore, we do not focus on studying the effect of horizons again from the perspective of algorithm, and the future horizon is set to be the same size as the past horizon (h
=
f) according to the previous conclusion. Then, there are three parameters (h
=
f, λ
1, and λ
2) to choose from. At first, only the effect of different horizons for the identification model is investigated, where other parameters (λ
1 and λ
2) remain unchanged. The truncation error, which equals the sum of the remaining retail of Markov parameters, is presented in Figure 5. And the identification result is shown in Figure 6.Figure 5
The truncation error against past horizonh.Figure 6
Effect of differenth and f (λ
1 and λ
2 are constant).As shown in Figures5 and 6, the identification precision is unsatisfactory when the horizon is small. In the case of a small horizon, the truncation error of Markov parameters is large; that means the approximation of A
-
h
=
0 is unsatisfied. And the outputs of the identified vehicle model are far from the test data. When the horizon is too big (h
>
40), the truncation error decreases slowly, but it would increase the calculation burden and show overshoot for the identified model. From the above results, it is found that the identified model is of high precision when h
=
f
=
40.The forgetting factor can overcome the data saturation in RLS method, and it is usually chosen between 0.9 and 1. In practice, the convergence rate is slow when a big forgetting factor is used, but the advantage is that the identified parameter is much more stable. When the forgetting factor is small, the algorithm is more sensitive to recent samples and the convergence rate is faster. However, the fast convergence rate is at the cost of sacrificing antijamming capability. Based on a vehicle pulse input test, the effect of forgetting factorsλ
1 and λ
2 is studied, respectively. And MSE (Mean Squared Error) listed in Table 2 is employed to evaluate the matching between test data and estimated outputs.Table 2
Results of MSEs under different forgetting factors.
λ
1
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
λ
2
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
0.92
0.3358
2.8668
0.92
0.9755
1.8686
0.95
1.4430
3.1832
0.95
0.2620
0.3311
0.98
2.9665
7.2588
0.97
0.5157
0.3658
0.992
0.2620
0.3311
0.99
0.2030
0.4169As observed in Figure7, when λ
1 is equal to 0.992, the identified model accords well with test data and the MSEs of both lateral acceleration and yaw rate are minimal. Considering that λ
1 is used to estimate the Markov parameter matrix Θ, a big λ
1 can also ensure the antijamming capability. As shown in Figure 8, the identified model is closer to the test data when λ
2 is 0.95 or 0.99. The MSE of lateral acceleration is smallest when λ
2
=
0.99. On the other hand, the MSE of yaw rate is best when λ
2
=
0.95. Balancing the convergence rate against high precision, a small forgetting factor can make the convergence rate increase. Therefore, λ
2
=
0.95 is a better choice.Figure 7
Effect ofλ
1 (h, f, and λ
2 are constant).Figure 8
Effect ofλ
2 (h, f, and λ
1 are constant).Through above analysis, a set of suitable parameters can be determined for the test vehicle. Whenh
=
f
=
40, λ
1
=
0.992, and λ
2
=
0.95, the identification method is of high precision. And this set of constant parameters is also used in the following test validation.


### 3.2. Vehicle Test Validation

In order to validate the vehicle identification model, the verification is carried out by a series of real vehicle tests. It is worth noting that the identified vehicle model at hand is a discrete-time model.Pylon course slalom tests at different vehicle speeds (30 km/h, 40 km/h, 50 km/h, and 60 km/h) are used to identify the vehicle system. As a contrast, the LTI vehicle model structure presented in (5) is identified by traditional subspace method, and the batch-wise identification algorithm in detail can be referred to in literature [9]. Due to the identified system matrices (A
(
t
), B
(
t
), C
(
t
), and K
(
t
)) being constantly updated at each step, it is not realistic to demonstrate all the elements of the matrices. A usual way adopted in [29, 31] is to show the eigenvalues of the identified matrix A
(
t
). The identified pole locations and the predicted outputs of identified models under 40 km/h and 60 km/h are given in Figures 9 and 10, respectively. Moreover, the MSEs listed in Table 3 can be used in quantitative analysis.Table 3
Results of MSEs under different vehicle speeds.
Speed(km/h)
LTV model
LTI model
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
30
0.009
0.162
0.021
0.280
40
0.016
0.176
0.074
0.806
50
0.029
0.309
0.280
1.456
60
0.061
0.446
0.817
2.822Figure 9
Identified pole location. “×”: the poles of LTV model; “+”: the poles of LTI model.
(a)
U = 40 km/h
(b)
U = 60 km/hFigure 10
Comparison between road test and different vehicle models under pylon course slalom test. (a) The average vehicle speed is 40 km/h. (b) The average vehicle speed is 60 km/h.
(a)
(b)As shown in Figure10, the LTI model can predict the vehicle behavior well when the average speed is 40 km/h. For the case of 60 km/h, the precision of the prediction outputs of LTI vehicle model decreases, where the MSEs of yaw rate are even over 2.8 (°/s)2. On the other hand, the predicted outputs of the LTV model always accord well with the real outputs, and the MSEs of lateral acceleration and yaw rate are 0.061 (m/s2)2and 0.446 (°/s)2, respectively, which is much smaller than that of LTI model. In fact, the precondition of the identifiable LTI model structure is the linear range for lateral acceleration, and this assumption is satisfied when a
y
<
0.4 g. With the increase of lateral acceleration, the time-varying characteristics of the vehicle intensify gradually. Therefore, it is necessary to investigate the effect of the proposed identification vehicle model under different lateral acceleration, especially big lateral acceleration.A series of step input tests with different lateral acceleration (0.2 g, 0.4 g, and 0.6 g) and a double lane change test with a maximum lateral acceleration of 0.8 g are exploited. The average vehicle speed is 100 km/h during these tests. As a contrast, the LTI vehicle model identified by traditional subspace method from step input data is also introduced. The predicted outputs of identified models are shown in Figures11 and 12. Moreover, the MSEs under different lateral acceleration are calculated in Table 4.Table 4
Results of MSEs under different lateral acceleration.
a
y (m/s2)
LTV model
LTI model
MSE ((m/s2)2)lateral acceleration
MSE ((°/s)2)yaw rate
MSE ((m/s2)2)lateral acceleration
MSE ((°/s)2)yaw rate
2
0.067
0.202
0.267
0.110
4
0.121
0.492
0.759
0.508
6
0.221
0.895
2.076
1.414
8
0.842
0.569
3.818
4.956Figure 11
Comparison among road test and different vehicle models under step input test.Figure 12
Comparison among road test and different vehicle models under double lane change test.As shown in Figures11 and 12, there is a little steady error for the predicted lateral acceleration of the LTI model in the cases of 0.2 g and 0.4 g, but the model is still in good agreement with the real outputs, and the acceptable MSEs of lateral acceleration and yaw rate are less than 0.075 (m/s2)2and 0.85 (°/s)2, respectively. Moreover, it is found that the LTI model becomes less effective in the case of 0.6 g and 0.8 g lateral acceleration. That is because the nonlinearity of actual vehicle system becomes prominent that the LTI model structure is no longer applicable, and the MSEs of LTI model become large. By contrast, the LTV model is satisfied. Even in the case of double lane change with the maximum lateral acceleration of 0.8 g, the LTV model can still predict the real vehicle outputs accurately, and the MSEs of lateral acceleration and yaw rate are only 0.842 (m/s2)2and 0.569 (°/s)2, respectively.


## 3.1. The Choice of Algorithm Parameters

Before vehicle identification, there are four constant parameters in the algorithm to be determined, and they are the past and future horizonsh and f and forgetting factors λ
1 and λ
2, respectively. These parameters play important roles in the convergence rate and identification precision. Based on vehicle pulse input test, the influences of the pending parameters are studied, and it could guide for an appropriate choice of these parameters.Generally, the past and future horizons are user choices in the subspace algorithm. And the future horizonf is less than or equal to the past horizon h (f
≤
h). As pointed out by Houtzager et al. [29], the truncation error decreased with the increase of past horizon h, but the computational time became larger. On the other hand, Chiuso [30] suggested that a large future horizon f was often better for identification when the signal-to-noise ratio was low. And Houtzager et al. [29] also indicated that when premultiplied weighting matrix was not applied for R
P
B
S
I
D
o
p
t algorithm, the high identification accuracy was obtained when the future horizon was close to the past horizon. Therefore, we do not focus on studying the effect of horizons again from the perspective of algorithm, and the future horizon is set to be the same size as the past horizon (h
=
f) according to the previous conclusion. Then, there are three parameters (h
=
f, λ
1, and λ
2) to choose from. At first, only the effect of different horizons for the identification model is investigated, where other parameters (λ
1 and λ
2) remain unchanged. The truncation error, which equals the sum of the remaining retail of Markov parameters, is presented in Figure 5. And the identification result is shown in Figure 6.Figure 5
The truncation error against past horizonh.Figure 6
Effect of differenth and f (λ
1 and λ
2 are constant).As shown in Figures5 and 6, the identification precision is unsatisfactory when the horizon is small. In the case of a small horizon, the truncation error of Markov parameters is large; that means the approximation of A
-
h
=
0 is unsatisfied. And the outputs of the identified vehicle model are far from the test data. When the horizon is too big (h
>
40), the truncation error decreases slowly, but it would increase the calculation burden and show overshoot for the identified model. From the above results, it is found that the identified model is of high precision when h
=
f
=
40.The forgetting factor can overcome the data saturation in RLS method, and it is usually chosen between 0.9 and 1. In practice, the convergence rate is slow when a big forgetting factor is used, but the advantage is that the identified parameter is much more stable. When the forgetting factor is small, the algorithm is more sensitive to recent samples and the convergence rate is faster. However, the fast convergence rate is at the cost of sacrificing antijamming capability. Based on a vehicle pulse input test, the effect of forgetting factorsλ
1 and λ
2 is studied, respectively. And MSE (Mean Squared Error) listed in Table 2 is employed to evaluate the matching between test data and estimated outputs.Table 2
Results of MSEs under different forgetting factors.
λ
1
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
λ
2
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
0.92
0.3358
2.8668
0.92
0.9755
1.8686
0.95
1.4430
3.1832
0.95
0.2620
0.3311
0.98
2.9665
7.2588
0.97
0.5157
0.3658
0.992
0.2620
0.3311
0.99
0.2030
0.4169As observed in Figure7, when λ
1 is equal to 0.992, the identified model accords well with test data and the MSEs of both lateral acceleration and yaw rate are minimal. Considering that λ
1 is used to estimate the Markov parameter matrix Θ, a big λ
1 can also ensure the antijamming capability. As shown in Figure 8, the identified model is closer to the test data when λ
2 is 0.95 or 0.99. The MSE of lateral acceleration is smallest when λ
2
=
0.99. On the other hand, the MSE of yaw rate is best when λ
2
=
0.95. Balancing the convergence rate against high precision, a small forgetting factor can make the convergence rate increase. Therefore, λ
2
=
0.95 is a better choice.Figure 7
Effect ofλ
1 (h, f, and λ
2 are constant).Figure 8
Effect ofλ
2 (h, f, and λ
1 are constant).Through above analysis, a set of suitable parameters can be determined for the test vehicle. Whenh
=
f
=
40, λ
1
=
0.992, and λ
2
=
0.95, the identification method is of high precision. And this set of constant parameters is also used in the following test validation.


## 3.2. Vehicle Test Validation

In order to validate the vehicle identification model, the verification is carried out by a series of real vehicle tests. It is worth noting that the identified vehicle model at hand is a discrete-time model.Pylon course slalom tests at different vehicle speeds (30 km/h, 40 km/h, 50 km/h, and 60 km/h) are used to identify the vehicle system. As a contrast, the LTI vehicle model structure presented in (5) is identified by traditional subspace method, and the batch-wise identification algorithm in detail can be referred to in literature [9]. Due to the identified system matrices (A
(
t
), B
(
t
), C
(
t
), and K
(
t
)) being constantly updated at each step, it is not realistic to demonstrate all the elements of the matrices. A usual way adopted in [29, 31] is to show the eigenvalues of the identified matrix A
(
t
). The identified pole locations and the predicted outputs of identified models under 40 km/h and 60 km/h are given in Figures 9 and 10, respectively. Moreover, the MSEs listed in Table 3 can be used in quantitative analysis.Table 3
Results of MSEs under different vehicle speeds.
Speed(km/h)
LTV model
LTI model
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
MSE ((m/s2)2)(lateral acceleration)
MSE ((°/s)2)(yaw rate)
30
0.009
0.162
0.021
0.280
40
0.016
0.176
0.074
0.806
50
0.029
0.309
0.280
1.456
60
0.061
0.446
0.817
2.822Figure 9
Identified pole location. “×”: the poles of LTV model; “+”: the poles of LTI model.
(a)
U = 40 km/h
(b)
U = 60 km/hFigure 10
Comparison between road test and different vehicle models under pylon course slalom test. (a) The average vehicle speed is 40 km/h. (b) The average vehicle speed is 60 km/h.
(a)
(b)As shown in Figure10, the LTI model can predict the vehicle behavior well when the average speed is 40 km/h. For the case of 60 km/h, the precision of the prediction outputs of LTI vehicle model decreases, where the MSEs of yaw rate are even over 2.8 (°/s)2. On the other hand, the predicted outputs of the LTV model always accord well with the real outputs, and the MSEs of lateral acceleration and yaw rate are 0.061 (m/s2)2and 0.446 (°/s)2, respectively, which is much smaller than that of LTI model. In fact, the precondition of the identifiable LTI model structure is the linear range for lateral acceleration, and this assumption is satisfied when a
y
<
0.4 g. With the increase of lateral acceleration, the time-varying characteristics of the vehicle intensify gradually. Therefore, it is necessary to investigate the effect of the proposed identification vehicle model under different lateral acceleration, especially big lateral acceleration.A series of step input tests with different lateral acceleration (0.2 g, 0.4 g, and 0.6 g) and a double lane change test with a maximum lateral acceleration of 0.8 g are exploited. The average vehicle speed is 100 km/h during these tests. As a contrast, the LTI vehicle model identified by traditional subspace method from step input data is also introduced. The predicted outputs of identified models are shown in Figures11 and 12. Moreover, the MSEs under different lateral acceleration are calculated in Table 4.Table 4
Results of MSEs under different lateral acceleration.
a
y (m/s2)
LTV model
LTI model
MSE ((m/s2)2)lateral acceleration
MSE ((°/s)2)yaw rate
MSE ((m/s2)2)lateral acceleration
MSE ((°/s)2)yaw rate
2
0.067
0.202
0.267
0.110
4
0.121
0.492
0.759
0.508
6
0.221
0.895
2.076
1.414
8
0.842
0.569
3.818
4.956Figure 11
Comparison among road test and different vehicle models under step input test.Figure 12
Comparison among road test and different vehicle models under double lane change test.As shown in Figures11 and 12, there is a little steady error for the predicted lateral acceleration of the LTI model in the cases of 0.2 g and 0.4 g, but the model is still in good agreement with the real outputs, and the acceptable MSEs of lateral acceleration and yaw rate are less than 0.075 (m/s2)2and 0.85 (°/s)2, respectively. Moreover, it is found that the LTI model becomes less effective in the case of 0.6 g and 0.8 g lateral acceleration. That is because the nonlinearity of actual vehicle system becomes prominent that the LTI model structure is no longer applicable, and the MSEs of LTI model become large. By contrast, the LTV model is satisfied. Even in the case of double lane change with the maximum lateral acceleration of 0.8 g, the LTV model can still predict the real vehicle outputs accurately, and the MSEs of lateral acceleration and yaw rate are only 0.842 (m/s2)2and 0.569 (°/s)2, respectively.


## 4. Discussions

Through above qualitative and quantitative analysis, the characteristics of the proposed vehicle identification model are discussed as follows.The advantages of the proposed LTV identification model mainly include the enhancement of practicability, improvement of the estimation precision, and the expansion of application scope. Firstly, compared with the existing LPV vehicle model structure, the proposed model is more practical. That is because the sideslip angle, an important output signal in LPV model structure, is costly obtained for most vehicles, while the proposed model structure is derived from the measurable signals, steering wheel angle, yaw rate, and lateral acceleration. Secondly, the estimation precision is improved in comparison with the LTI vehicle model identified by traditional subspace method. Validated by a series of standard road tests, the identified LTV model predicts the vehicle outputs more accurately, and the MSEs of LTV model are satisfactory under different maneuvers. At last, the application scope of the proposed identification model is not restricted for small lateral acceleration. Due to the vehicle nonlinearity becoming severe and the LTI vehicle model being inappropriate for the limitation of model structure when lateral acceleration is over 0.4 g, the effect of the LTI identified model gets worse. By experimental verification, the proposed identification model is proved applicable and accurate enough to predict vehicle behavior even when the maximum lateral acceleration is 0.8 g.


## 5. Conclusions

In the research of vehicle modeling by identification, a 4th-order identifiable LTV vehicle model structure is proposed to describe real vehicle lateral dynamics behavior. To identify the LTV vehicle model, a recursive subspace identification method combining theP
B
S
I
D
o
p
t and OPAST technology is introduced, which can provide asymptotically consistent estimation and have a global convergence property. Based on pulse input test data, the effects of four parameters in the algorithm are studied, and a set of practical optimal parameters is selected for the vehicle model identification. By using a series of real vehicle tests under different maneuvers, the linear time-varying vehicle model can be identified in real-time. In comparison with identified LTI vehicle model, the proposed vehicle model can predict the real vehicle behavior accurately. Even when the lateral acceleration is over 0.4 g, the severe nonlinearity in vehicle system, the LTV vehicle model is still of high precision.Summarizing the paper, noteworthy contributions of this work are primarily twofold: an identifiable LTV vehicle model structure which is an improvement to the existing model structure and the development of a practical subspace identification method to vehicle modeling. These contributions are validated by qualitative and quantitative analysis based on real vehicle tests. In fact, the proposed vehicle model is not suitable for intense cases, such as rapidly changing speed and fast steering. And the derivation calculation for steering angular velocity as an input signal is not convenient. In future work, a more general model to describe the vehicle system under extreme working conditions should be studied.



---

*Source: 1715762-2016-06-22.xml*