Application of Recursive Subspace Method in Vehicle Lateral Dynamics Model Identification

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.


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   or 1/V  .To identify the multiinput multioutput (MIMO) vehicle system, Cerone et al. [8] addressed a set-membership method, but the complex and low efficiency identification process consisted of two singleinput 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 PBSID opt , 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 PBSID opt 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.

Vehicle Model Identification Method
where  is the vehicle mass, V  is the longitudinal speed, V  is the lateral speed of the CG of the vehicle,  is the yaw rate,   is the vehicle yaw moment inertia, and   and   are front and rear tire lateral forces (Figure 1).
For linear tire model under small tire slip angle, the lateral forces can be written as where   and   are the average cornering stiffness of front and rear axles and   and   are the front and rear tire slip angles, respectively.The tire slip angles can be expressed approximatively as Substituting ( 2) and ( 3) into (1), the vehicle state equation can be described as In practice, only the lateral acceleration   and yaw rate  are available by vehicle sensors, but the measurable outputs cannot be expressed by the linear combination of state variables V  and .In fact, the measurable output   is equal to the sum of vehicle speed V  multiplied by yaw rate  and the derivative of lateral velocity V  (  = V   + V  ).Therefore, it is necessary to extend (4) so as to obtain the state variable V  .Taking the derivative of (4), an extended state equation, including the derivative of lateral velocity, can be obtained as follows: Given a longitudinal speed V  , 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.
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: where , and system matrices are ) , 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 the RPBSID opt method is that the system to be identified should be minimal, and this means system (6) should be observable and controllable for any time .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   (),   (), and V  () 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 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.

Review of 𝑅𝑃𝐵𝑆𝐼𝐷 𝑜𝑝𝑡 Algorithm.
The PBSID opt method introduced in [13,22] is becoming an effective and wellestablished 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 PBSID opt (RPBSID opt ) identification.
Considering the process and measurement noises, the discrete form of an LTI system given can be expressed as where u() ∈ R  and () ∈ R  are input and output sequences, respectively.x() ∈ R  is a state vector, and e() ∈ R  is defined as e() = y() − ŷ( |  − 1), where ŷ( |  − 1) means the estimation of y() by using the data before moment −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 dimension  of x() is as small as possible.
(A2) The matrix A − KC is asymptotically stable.
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 ℎ 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 ℎ 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 matrices A, B, C, K, and the innovation process e() by the measurements of input and output (I/O) sequences {u(), y()} = =1 .The system expressed by (10) and (11) can be rewritten in a predictor form as process of the RPBSID opt 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 According to the assumption (A2), A ℎ approaches to 0 when the past horizon ℎ 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: where Iteration of ( 12) and ( 13) results in the subspace matrix formulation of the following form as where ] , The extended observability matrix is given as )  ]  , and the lower triangular block-Toeplitz matrix H  is given by According to (15), the state matrix X  can be expressed as where Substituting ( 19) into ( 16), then The matrix M is defined by where m  := CA −1 B ( = 1, 2, . ..) is called the predictor Markov parameters.As A ℎ = 0 whenever ℎ is sufficiently Mathematical Problems in Engineering large, m  = 0 ( > ℎ).Then, the matrix M in ( 21) can be simplified to the following form: According to ( 16) and ( 20), an important identity for PBSID opt algorithm can be found as below: And the first column vector form of ( 23) is The identity is the key equation for PBSID opt algorithm.
Step 1 (the estimation of Markov parameters m  ).Substitution of ( 15) into ( 13) results in the following form: According to (13), the output vector y() can be expressed as where , and the element m  := CA −1 B ( = 1, . . ., ℎ) is called predictor Markov parameters: The Markov parameters m  can be estimated by RLS (Recursive Least Square) filter with forgetting factor  1 .
Step 2 (the update of the extended observability matrix Γ  ).An OPAST method is proposed to update the extended matrix Γ  .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  as follows: where  2 is a forgetting factor and W ∈ R × ( > ) is a full rank matrix.
Theorem 1 (see [19]).The global minimum (W) is attained if and only if W = QT, where Q contains the  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 if R  is replaced with The key idea of PAST is to replace W  ()() in (27) with The projection approximation is valid if the weighted matrix is slowly varying.Then the criterion of quadratic form with respect to W() is rewritten as follows: which is minimized by where Rℎ () = ∑  =1  − 2 ()h  () and Rℎ () = ∑  =1  − 2 h()h  () are covariance matrices, and it is assumed that the inverse of matrix Rℎ () exists.
Since the criterion (W) is used for approximation of (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: As observed in identity equation (24), the column space of Mz ℎ () is consistent with the column space of Γ  .Mz ℎ () is fed into objective function (27) as the vector , and then the output W() is corresponding to the extended observability matrix Γ  .And this rationale is widely used in previous literature [17,28] for updating observability matrix Γ  .
Step 3 (the estimation of the system matrices {A, B, C, K}).Calculate the estimated matrix Ĉ.It can be extracted directly and equals the front  rows of Γ  according to the definition: where Γ  (1 : , :) is a MATLAB expression, and it indicates the front  rows.

Input vector
Output vector

Real vehicle system
Step 1 Step 2 Step 3

Identification module
Predicted outputs Identified vehicle model matrixes where e(k) = y(k a y (t) According to (10), the system matrices { Â() B() K()} can be estimated by minimizing the following criterion: where the innovation e() = y() − Ĉ()x().And the states x can be estimated according to (24), where matrices Γ  and M are obtained before and the operation " †" donates Moore-Penrose inverse.

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 horizons ℎ and  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 horizon  is less than or equal to the past horizon ℎ ( ≤ ℎ).As pointed out by Houtzager et al. [29], the truncation error decreased with the increase of past horizon ℎ, but the computational time became larger.On the other hand, Chiuso [30] suggested that a large future horizon  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 RPBSID opt 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 (ℎ = ) according to the previous conclusion.Then, there are three parameters (ℎ = ,  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.
As shown in Figures 5 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 ℎ = 0 is unsatisfied.And the outputs of the identified vehicle model are far from the test data.When the horizon is too big (ℎ > 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 ℎ =  = 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.
As observed in Figure 7, 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.Through above analysis, a set of suitable parameters can be determined for the test vehicle.When ℎ =  = 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.

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(), B(), C(), and K()) 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().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.
As shown in Figure 10, 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/s 2 ) 2 and 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   < 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 Figures 11 and 12.Moreover, the MSEs under different lateral acceleration are calculated in Table 4.
As shown in Figures 11 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/s 2 ) 2 and 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/s 2 ) 2 and 0.569 ( ∘ /s) 2 , respectively.

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

Conclusions
In the research of vehicle modeling by identification, a 4thorder 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 the PBSID opt 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.

Appendix
According to (8) In consideration of safe driving, most of real vehicles are understeer, which means   −   > 0. It can be seen that the elements in the 3rd and 4th rows are not equal to 0, and the first 4 rows of matrix Γ are linearly independent (rank(Γ) = 4).Thus, the system is observable.According to (9), the controllability matrix Δ is

Figure 2 :
Figure 2: The relationship between tire lateral force and slip angle (  is the vertical load).

Figure 4 :
Figure 4: The vehicle and instruments in road test.

Figure 10 :
Figure 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.

Figure 12 :
Figure 12: Comparison among road test and different vehicle models under double lane change test.

Table 1 :
List of the characteristics of instruments in the road test.one of the output signals.The lateral acceleration   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 (6)cess.Due to the timevarying 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 =  G and B = ∫  0  G H , where  is a sampling interval.Through above approximation, the RPBSID opt method can be applied to estimate and update the model matrices at each moment , and the detailed process is given in Figure3.As 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 sensor cluster of accelerometer and gyroscope.And a lowpass Kaiser Filter is used for the collected signals.During the road tests, sample frequency is set as 100 Hz.And Table1lists the characteristic of the instruments in the road tests.

Table 2 :
Results of MSEs under different forgetting factors.

Table 3 :
Results of MSEs under different vehicle speeds.

Table 4 :
Results of MSEs under different lateral acceleration.