Short Data-Based Output-Only Identification for Time-Varying Systems with Fast Dynamic Evolution

,


Introduction
e need for time-varying (TV) system characterization is pervasive throughout the aerospace, mechanical, transportation, and manufacturing engineering practice.It is usually difficult to build the explicit dynamic model of a practical system via mechanism analysis; hence, TV system identification in a way that takes time variation explicitly into account is receiving more and more attention [1][2][3][4].With the development of engineering applications, systems with fast varying dynamics are being widely used and their intrinsic nonstationary characteristics are increasingly inevitable.On the one hand, these systems may change appreciably over a relatively short time interval due to their fast evolution in the dynamics.On the other hand, the dynamic characteristics of TV systems may not fit into a laboratory and controlled testing may not be feasible under operating conditions.In other words, TV system's ambient excitation is usually difficult to measure under operating conditions and its dynamics have to be determined without measuring the excitation.erefore, short data-based output-only identification for TV systems is considered in this work in order to capture their fast varying dynamics.
Identification methods for TV systems are generally classified under the umbrella of time-frequency methods [1].Most frequency-domain methods employ time-frequency analysis (e.g., the short time Fourier transform, Cohen's class, wavelet-based methods, the Hilbert-Huang transform, etc.), which are based on nonparametric representations of the nonstationary signal as a simultaneous function of time and frequency.By contrast, time-domain methods employ parametric representations (e.g., time-dependent autoregressive moving average (TARMA) models), which offer a number of potential advantages, such as improved accuracy, resolution, and tracking of the TV dynamics.TARMA model-based methods may be further divided into three classes according to the type of "structure" imposed upon the evolution of the TV model parameters: unstructured parameter evolution (UPE), stochastic parameter evolution (SPE), and deterministic parameter evolution (DPE) methods [3]. e DPE methods are mainly based on explicit models of parameter evolution through projecting the parameter trajectory onto specific functional subspaces defined by deterministic basis functions (such as trigonometric, Legendre, Chebyshev, wavelets, B-splines, moving Kriging shape functions, radial basis functions, and others).By postulating TV model parameters as deterministic functions of time, the TARMA model can be transformed into a functional series (FS) TARMA model with time-invariant parameter matrix.erefore, DPE methods are known to track fast dynamic evolution by selecting proper functional subspaces over their UPE and SPE counterparts [2][3][4].
DPE methods have played an important role in the development of TV system identification and nonstationary signal processing.In 1970, Rao [5] approximated AR coefficients by truncated Taylor series expansion and introduced the idea of representing TV model parameters as polynomials of time.ereafter, Mendel [6] summarized the nonstationary identification methods and classified them into UPE, SPE, and DPE, which are still in use today.Kozin [7] approximated AR coefficients by a linear combination of Legendre basis functions and applied the FS-TAR model to earthquake ground motion analysis.In 1983, Grenier [8] extended the FS-TAR model to the FS-TARMA case by examining the existence and uniqueness of a FS-TARMA model representation for a nonstationary system, which was successfully applied to speech signal modeling [9][10][11][12].Gersh and Kitagawa [13] extended the univariate FS-TAR model to the multivariate/vector case.Niedzwiecki [14,15] defined the idea of expanding TV model parameters onto linear combinations of basis functions as "functional series modeling" and proposed a recursive least squares for the estimation of FS-TAR model parameters.Tsatsanis and Giannakis [16,17] proposed FS-TARMA models with exogenous input and FS time-dependent finite impulse response models and considered the issue of model structure selection.In 2000, Niedzwiecki [2] summarized the DPE methods and investigated the time and frequency characteristics of basis function estimators.
In 2006, Poulimenos and Fassois [3] systematically reviewed the UPE, SPE, and DPE methods in terms of achievable model parsimony, model parameter estimation accuracy, tracking capability of TV dynamics, and computational simplicity.Spiridonakos and Fassois extended the FS-TAR model to the FS-TARMA case [18] and proposed the adaptable FS-TARMA model characterized by basis functions that are adaptable to the signal being modeled [19].In 2014, Spiridonakos and Fassois [4] reviewed the conventional and adaptable FS-TARMA model in terms of parameter estimation, model structure selection, and model validation and outlined the developments of the DPE methods.Recently, Yang et al. [20] have proposed an adaptable FS-TARMA modeling method by selecting moving Kriging shape functions as basis functions.Zhou et al. [21] have proposed a least squares support vector machine-based FS-TAR modeling method by selecting compactly supported radial basis functions.Bertha and Golinval [22] have focused on the identification of TV mode shapes during the FS-TARMA modeling.Li et al. [23] have focused on the problem of extracting the real modal parameters from computational ones during the adaptable FS-TAR modeling.Ma et al. [24,25] have focused on the recursive identification of TV systems by introducing the idea of kernel methods into the FS-TARMA modeling.
Despite the foregoing progress that has been achieved so far, there are still many issues that are open and require work, for example, the feasibility of TV system identification under a relatively limited data length.DPE methods can achieve better model parsimony over their UPE and SPE counterparts, yet they still require a significant number of parameters to represent complicated time-dependencies and dynamics characterized by larger numbers of degrees-offreedom [3].erefore, the aim of the paper is two-fold: (1) the existing least squares-based FS-TAR modeling is reviewed to reveal the problem of TV system identification under a relatively limited data length and (2) a ridge regression-based FS-TAR method is subsequently proposed to track fast dynamic evolution of the system being modeled by using short data.
e paper is organized as follows: the least squares-based FS-TAR modeling is presented in Section 2. Section 3 proposes a new ridge regression-based FS-TAR method to identify TV systems with fast dynamic evolution under limited data length.e proposed and existing methods are numerically and experimentally tested in Section 4 and Section 5, respectively.Finally, Section 6 gives the remarkable conclusions.

LS-Based FS-TAR Model
Parameter Estimation

FS-TAR Model.
A time-dependent autoregressive (TAR) model with n a designating its autoregressive (AR) order is given by [18] x Any given function can be approximated by a set of selected basis functions with arbitrary accuracy, as long as a sufficient number of basis functions is used [26].erefore, a TV parameter v[t] can be approximated by 2 Shock and Vibration where G b v (j) [t] denotes a set of basis functions selected from a suitable functional subspace, v j the coefficient of projection, p v the subspace dimensionality, and the index b v (i), i � 1, . . .p v the specific basis functions of the selected functional subspace that are included in the basis.erefore, the time-dependent AR parameter matrix A i [t] can be expressed as [3,4,18,20] e TAR model can be rewritten into the following FS-TAR model by substituting equation (3) into equation ( 1): being the time-invariant parameter matrix, and the corresponding regression vector, where ⊗ is the Kronecker product [27], g b a [t] is given by Evidently, the time-varying TAR model of equation ( 1) is transformed into a time-invariant FS-TAR model by expanding the time-dependent model parameters onto the selected functional subspaces.As given by equation ( 5), the FS-TAR model is fully described by k 2 n a p a parameters (coefficients of projection).

Least Squares Estimation.
Suppose the length of measured response data is N, the matrix equation can be obtained by assembling each equation of form (4) for t � 1, 2, . . ., N, as follows: with , and e parameter matrix θ can be estimated by using the well-known least-squares estimator through minimizing the sum of squared residuals, as follows: where (•) + denotes the Moore-Penrose pseudo inverse.e time-dependent AR parameter matrix A i [t] can be subsequently estimated by using equation ( 3), after substituting θ by the obtained estimate  θ.However, if measured response data length is short and the inequality N < k 2 n a p a holds, there are fewer equations than unknown parameters and more than one θ satisfies equation (8).In other words, the solution is not unique and the problem is ill posed.In such cases, the ordinary least squares estimation is sensitive to errors and may lead to unreliable estimates.From the standpoint of the FS-TAR modeling, the data can be considered as "short" if the sample length (N) is smaller than the number of unknown parameters of the identified model (k 2 n a p a ).

RR-Based FS-TAR Model Parameter Estimation
3.1.Ridge Regression Estimation.Ridge regression, also known as Tikhonov regularization, is a commonly used method of regularization of ill-posed problems [28,29].
Given the regularization parameter c ≥ 0, the least squares cost function is defined as By using the ridge regression, we have the solution to equation (10) as where By using the matrix inversion lemma, we have Shock and Vibration Obviously, the computational complexity of the matrix inversion in equations ( 11) and ( 12) is, respectively, O(k 3 n 3 a p 3 a ) and O(t 3 ).In order to avoid the matrix inversion, we further define Q[t] and k[t] as en, we have e update rule for the inverse of the above growing matrix is By combining equations ( 12), (13), and ( 16), we have with the final estimate  θ �  θ[N].Similarly, the time-dependent AR parameter matrix A i [t] can be subsequently estimated by using equation (3), after substituting θ by the obtained estimate  θ.

FS-TAR Model-Based Structural Dynamics. Once A i [t]
have been estimated, the system's "frozen-time" poles can be obtained by solving the following general eigenvalue problem [20]: where I kn a designates the identity matrix of order kn a and p r [t] and T , respectively, the rth poles and eigenvector of D[t] at time instant t with the rth mode shape L r .e matrix D[t] is constructed from e "frozen-time" natural frequencies of the system can be computed by where | • | designates the absolute value and Δt the sampling period.
e response "frozen-time" power spectral density (PSD) for each time instant can be obtained as where ω designates frequency in rad/s, j the imaginary unit, (•) H the Hermitian transpose, and A[e −jωΔt , t] �

Numerical Validation
e performance characteristics of the proposed RR-FS-TAR method, along with the LS-FS-TAR method, are examined via Monte Carlo study focused on the identification of a TV three degrees-of-freedom structural system subject to a Gaussian, zero-mean, and uncorrelated excitation r(t), as shown in Figure 1 [3,30].Numerical values of the system parameters are given in Table 1.
ree cases are considered by choosing different values of the scaling factor N s .Case 1: slow dynamic evolution with N s � 0.2, Case 2: medium dynamic evolution with N s � 1, and Case 3: fast dynamic evolution with N s � 5. e simulation time is, respectively, 425, 85, and 17 seconds for each case.Given a particular excitation realization, the displacement responses are obtained through the Runge-Kutta 4/5 method characterized by variable integration step and subsequently recorded at a sampling frequency 12 Hz.erefore, the length of the samples is, respectively, 5100, 1020, and 204 for each case.e above procedure is repeated for each Monte Carlo test (number of tests: 50), with a different excitation realization used each time.

Case 1: Slow Dynamic Evolution.
In this case, the displacement responses are 425s long and the sample length is N � 5100.
e AR order is manually selected as n a � 2 because the actual number of modes present is N m � 3 and n a � 2N m /k � 2, as given by equation (20).e functional basis spanned by trigonometric functions of the form [3,30], is selected.e subspace dimensionality p a � 40 and index b a � [0 : 39] are selected to approximate the TV model parameters.Figure 2 depicts the LS-based and RR-based natural frequency estimates from the 50 Monte Carlo tests and their true counterparts.As shown in Figure 2 the TV natural frequencies are adequately tracked by the proposed RR-FS-TAR method and the LS-FS-TAR method, due to the relatively larger sample length.However, the proposed RR-FS-TAR method seems to attain better overall performance in terms of tracking and accuracy as the LS-FS-TAR method- Table 1: Numerical values of the system parameters.Shock and Vibration quite rough.e trend of the natural frequency estimates obtained by the LS-FS-TAR method is in overall agreement with the true values, but they are of rather poor quality in terms of tracking and accuracy.By contrast, the proposed RR-FS-TAR method attains much better overall performance in terms of tracking and accuracy, even if the system changes faster and the available sample length is relatively shorter.

Case 3: Fast Dynamic Evolution.
In this case, the system changes appreciably over a relatively short time interval.e displacement responses are 17s long and the sample length is N � 204.e AR order is selected as n a � 2. e functional basis spanned by trigonometric functions of equation ( 22) is here selected with the same values of subspace dimensionality p a � 40 and index b a � [0 : 39]. Figure 4 depicts the LS-based and RR-based natural frequency estimates from the 50 Monte Carlo tests and their true counterparts.e sample length (N � 204) is much smaller than the number of unknown parameters (k 2 n a p a � 720).In such case, the ordinary least squares estimation is sensitive to errors and may lead to unreliable estimates.is behavior is confirmed by Figure 4, where the TV natural frequencies are not tracked by the LS-FS-TAR method and their estimates are rather vaguely shown or even missing.e proposed RR-FS-TAR method is able to deal with the ill-posed problem via regularization.As shown in Figure 4, the natural frequency estimates obtained by the RR-FS-TAR method is in overall agreement with the true values, although they are rough and exhibit larger tracking errors.

Power Spectral Density Estimates.
e superior performance of the proposed RR-FS-TAR method is pronounced by the comparison of Figures 2-4 where the estimated natural frequencies are shown.In this section, we further examine the enhanced capability of the RR-FS-TAR method in terms of achievable power spectral density tracking and resolution.Figure 5 depicts the LS-based and RR-based power spectral density estimates from the single test for the TV numerical system with slow, medium, and fast dynamic evolution.
e LS-based and RR-based power spectral density estimates are quite similar for the TV numerical system with slow dynamic evolution, as shown in Figure 5.However, with the decrease of the sample length, the LS-FS-TAR method obtains worse results, while the proposed RR-FS-TAR method attains much better performance in terms of tracking and resolution, even if the system changes very fast and the available sample length is quite short (N � 204).In contrast to LS-based FS-TAR estimator requiring larger sample length to represent complicated time-dependencies and dynamics, RR-based FS-TAR estimator is able to track fast dynamic evolution of the system being modeled by using short data.

Experimental System.
e experimental system [25,31] is composed of the TV structure, the exciter system (an exciter and a power amplifier), force and motion transducers (an impedance head and 15 accelerometers), the measurement and acquisition module (LMS SCADAS III system), control systems (two Faulhaber DC motors and corresponding motion controllers), and boundary conditions, as shown in Figure 6. e TV structure is a coupled system consisting of a simply supported beam and a moving mass.e mass, driven by two DC motors, can slide on the beam.In this section, two cases are considered by adjusting the moving speed of the mass, as follows: Case 1: the speed is 0.130 m/s and Case 2: the speed is 0.325 m/s.
During the test, the mass starts at 0.5 m and stops at 1.8 m. erefore, the durations for the mass to move over 1.3 m are, respectively, 10 and 4 seconds.e shaker is switched off and removed, and the TV structure is only excited by the motion of the moving mass.Obviously, the interaction between the moving mass and the beam is  erefore, the excitation of the experimental system is unobservable and the system's dynamic characteristics have to be determined by exclusively using measured vibration responses.e vibration responses of the simply supported beam are measured by accelerometers, placed at 15 uniformly distributed positions along the beam from left to right, as shown in Figure 6. e acceleration responses are measured at a frequency 512 Hz and subsequently resampled at 256 Hz. erefore, the length of the samples is, respectively, 2560 and 1024 for the two cases.e above procedure is repeated for each Monte Carlo test (number of tests: 30), and the TV structure undergoes the same variation in each test.

Identification Results and Discussions.
e proposed RR-FS-TAR method and the LS-FS-TAR method are experimentally compared and assessed based on the nonstationary acceleration measurements of the laboratory TV structure.
e AR order is manually selected as n a � 2, and the functional basis spanned by trigonometric functions of equation ( 22) is selected with the subspace dimensionality p a � 20 and index b a � [0 : 19].Obviously, the sample length (N � 2560 or 1024) is much smaller than the number of unknown parameters (k 2 n a p a � 9000) and the least squares estimation may lead to unreliable estimates.
Figure 7 depicts the LS-based and RR-based natural frequency estimates from the 30 Monte Carlo tests when N � 2560.e first four baseline natural frequencies of the TV structure are obtained by using the "frozen-time" technique [31] and denoted by green dashed lines in Figure 7. e gray dots denote the natural frequency estimates from the 30 tests, and the red dots denote the natural frequency estimates from the single test.Evidently, the TV natural frequencies are not adequately tracked by the LS-FS-TAR method and their estimates are quite rough, while the proposed RR-FS-TAR method attains much better tracking performance and the corresponding natural frequency estimates exhibit less false modes.e LS-based and RR-based power spectral density estimates from the single test for the TV structure are also shown in Figure 7. Compared to the LS-based power spectral density estimates, the RR-based estimates appear smoother and are of much better quality in terms of tracking and resolution.
Similarly, Figure 8 depicts the LS-based and RR-based natural frequency estimates from the 30 Monte Carlo tests and power spectral density estimates from the single test when N � 1024.On the one hand, with the decrease of the sample length from 2560 to 1024, both the LS-FS-TAR and RR-FS-TAR methods obtain worse results.Compared to the identification results in Figure 7, natural frequency estimates exhibit larger overall scatter and power spectral density estimates appear more rough, as shown in Figure 8.On the other hand, the enhanced capability of the proposed RR-FS-TAR method is further pronounced to track fast dynamic evolution of the system based on short data.Compared to the LS-based identification results, the RR-based natural frequency estimates exhibit lower overall scatter, less false modes, and smaller tracking errors.Furthermore, the RRbased power spectral density estimates are quite smooth and the ridges are in good agreement with the baseline values, even if the structure changes very fast and the available sample length is relatively short.
For practical systems with larger numbers of degrees-offreedom and complicated TV dynamics, the number of unknown parameters (k 2 n a p a ) of the FS-TAR model is usually quite large.In the meantime, the measured response data length (N) is limited, when systems change appreciably over a relatively short time interval.In such cases, the ordinary least squares estimation may lead to less accurate or even unreliable estimates, and the proposed ridge regression-based FS-TAR method is suggested to identify these systems with fast dynamic evolution based on short data.Shock and Vibration

Conclusions
e FS-TAR model can achieve good representation parsimony, yet it requires a signi cant number of parameters to represent complicated time-dependencies and dynamics characterized by larger numbers of degrees-of-freedom.erefore, relatively longer data are usually needed to overcome ill-posed problems and ensure the reliability of least square estimation.However, in many practical cases, systems change appreciably over a relatively short time interval and their dynamics have to be determined based on limited measurement data.In this work, a short data-based output-only identi cation method for TV systems with fast dynamic evolution was proposed via regularization of illposed problems.e proposed ridge regression-based FS-TAR method was employed to identify numerical and experimental systems with di erent dynamic evolution and compared against the least squares-based FS-TAR method via Monte Carlo study.Comparison results con rm the advantages of the proposed method in terms of achievable natural frequency and power spectral density tracking, accuracy, and resolution of TV systems with fast dynamic evolution, when the available measurement data length is relatively short.
Data Availability e data used to support the ndings of this study are available from the authors upon request.

4
Shock and Vibrationbased natural frequency estimates exhibit some local tracking errors.4.2.Case 2: Medium Dynamic Evolution.In this case, the displacement responses are 85s long and the sample length is N 1020.e AR order is selected as n a 2. e functional basis spanned by trigonometric functions of equation (22) is here selected with the same values of subspace dimensionality p a 40 and index b a [0 : 39].Figure3depicts the LS-based and RR-based natural frequency estimates from the 50 Monte Carlo tests and their true counterparts.Evidently, the TV natural frequencies are not adequately tracked by the LS-FS-TAR method and their estimates are

Figure 2 :
Figure 2: Natural frequency estimates by LS (a) and RR (b) for the TV numerical system with slow dynamic evolution (gray dots: natural frequency estimates from the 50 Monte Carlo tests; red dots: natural frequency estimates from the single test; green dashed lines: true values of natural frequencies).

Figure 3 :
Figure 3: Natural frequency estimates by LS (a) and RR (b) for the TV numerical system with medium dynamic evolution (gray dots: natural frequency estimates from the 50 Monte Carlo tests; red dots: natural frequency estimates from the single test; green dashed lines: true values of natural frequencies).

6
Shock and Vibration difficult to measure.

Figure 4 :
Figure 4: Natural frequency estimates by LS (a) and RR (b) for the TV numerical system with fast dynamic evolution (gray dots: natural frequency estimates from the 50 Monte Carlo tests; red dots: natural frequency estimates from the single test; green dashed lines: true values of natural frequencies).

Figure 6 :Figure 5 :
Figure 6: Schematic diagram of the experimental system and photo of its setup.

Figure 7 :
Figure 7: Natural frequency and power spectral density estimates by LS (a, c) and RR (b, d) for the laboratory TV structure when N 2560 (gray dots: natural frequency estimates from the 30 Monte Carlo tests; red dots: natural frequency estimates from the single test; green dashed lines: baseline values of natural frequencies).

Figure 8 :
Figure 8: Natural frequency and power spectral density estimates by LS (a, c) and RR (b, d) for the laboratory TV structure when N 1024 (gray dots: natural frequency estimates from the 30 Monte Carlo tests; red dots: natural frequency estimates from the single test; green dashed lines: baseline values of natural frequencies).