Identification of linear time-varying dynamic systems based on the WKB method

This work proposes a new parameter identification method based on the Wentzel-Kramers-Brillouin (WKB) approximation for slow linear time-varying (LTV) dynamic systems. The considered time period is divided into a series of short time windows. In each time window, the assumption of “short time linearly varying” parameters is employed, and a nonlinear optimization problem is solved using the WKB results for the slow LTV dynamic system. A search algorithm is developed to find the optimal solution. In the identification process, only one type of response signal (displacement, velocity or acceleration) is required. Thus, numerical differentiation or integration of the measured signal, which leads to truncation or cumulative errors in noise environment, is avoided. The accuracy and robustness of the new identification method are validated by applying it to a particular LTV system with time-varying stiffness.

numerical differentiation and integration of one type of the signals that is measured, existing noise can lead to remarkable errors due to the truncation or cumulative errors [19]. To improve existing methods fundamentally, it is essential to introduce an analytical expression of response solution to LTV system identification.
One common way to solve identification problem for a LTV system is to use the state space model in discrete-time form [20][21][22][23]. However, the displacement and velocity data is essential in the recursive process, which limits application of this approach to the inverse problems when these signals cannot be obtained simultaneously. From this perspective, the Wentzel-Kramers-Brillouin (WKB) approximation is a promising tool for LTV system identification. The WKB approximation was proposed in 1920s to solve the 1D stationary Schrödinger equation in quantum mechanics with explicit approximate analytical representations [24]. Recently, it has been applied to problems on deformable solids, such as the analysis of transverse free vibrations of variable-section beams [25], catenary-shaped slender structures [26] and other [27][28][29][30], which can be reduced to solving ordinary differential equations (ODE) with variable coefficients.
In this research, the WKB method is adopted to obtain an approximate response for the linear system with slow TV stiffness under external excitation. In particular, given the external excitation, a specific functional relation between the displacement response and the TV stiffness of the system can be obtained using the WKB method. Subsequently, relations between the velocity or acceleration of the system and the TV stiffness can also be derived from the displacement function. Besides, to effectively track the change of the stiffness, the "short time linearly varying" assumption which was proposed by Chen et al. to identify a nonlinear TV system is adopted in this work [19]. By describing a TV parameter with linear temporal functions in different short time windows, this assumption has been proved to have a high identification efficiency and good robustness to the measurement noise [31].
Then, combined with the "short time linearly varying" assumption, for any type of the response signals, a nonlinear optimization problem can be formulated, in which the linearly varying stiffness minimizing the mean absolute errors is to be determined. To find an optimal solution for the TV stiffness, a search algorithm is developed. A numerical example is considered in this paper to validate the proposed identification method. The TV stiffness is identified for cases with different levels of measurement noise, with each signal from displacement, velocity and acceleration. To show the superiority in accuracy and robustness, the traditional identification method based on linear regression equations that requires all the response signals is also employed for comparison.
The main novelty of this research lies in the introduction of the WKB method to the LTV system identification. The shortcoming of the linear regression problem in Ref. [19] can be overcome by the explicit form of the WKB solution. Besides, to enable the use of the WKB method for LTV system identification, specific modifications are made to the conventional WKB solution, and an efficient search strategy is proposed for solving the parameters in the nonlinear optimization.
The rest of the paper is organized as follows. In Sect. 2, the WKB approximate solution for a linear dynamic system with TV stiffness is briefly derived. In Sect. 3, two critical issues in LTV system identification are discussed. The new LTV system identification method using a modified WKB method combined with "short time linearly varying" assumption is proposed in Sect. 4. To estimate the unknown TV stiffness accurately and efficiently, a new nonlinear search algorithm is proposed for the optimization problem. Besides, the identification process using velocity or acceleration signal is also considered in this section. Some numerical results of the implementation of the new identification method are discussed in Sect. 5. Finally, Sect. 6 gives the conclusive remarks and future directions.

WKB method for response prediction of LTV dynamic systems
In this section, the basics of the WKB method for LTV dynamic systems [19] is introduced. For the convenience of explaining the idea and principle of this work, we focus on the single degree of freedom (SDOF) system. The theoretical derivation is based on the SDOF system in this paper. The considered governing equation of motion of a SDOF system with a TV stiffness can be written as where c, k(t), and f are damping, TV stiffness and the external force respectively; x is the displacement and '.' and '..' denote the first and second derivatives with respect to time t. Equation (1) cannot be solved directly analytically. Numerical methods for ODEs, such as the Runge-Kutta method, are generally used to calculate the system response. Instead, in the present research, the WKB method is employed to obtain an approximate solution of Eq. (1) in an explicit form. As mentioned in [32], a system is slowly varying if its coefficients variation is much slower than the solution variation. Thus, the quantitative measure for the slowness ratio used in this paper is wherek is the time-derivative of k, is the frequency of external excitation and λ is the natural frequency of the system. Under the assumption of slowly varying parameters, Eq. (1) can be rewritten as in which the stiffness k varies in a slow time scale τ εt, where ε is a small parameter (0 < ε 1). Then, Eq. (3) is rewritten with respect to the independent time variable τ with d dt ε · d dτ . We assume that the LTV system satisfies the following conditions at the start τ s and the end τ e of the time interval [τ s , τ e ] The WKB approximation of a dynamic system with TV mass, damping and stiffness has been given in Ref. [33]. As a special case containing a TV stiffness, the theoretical derivation of the WKB solution of Eq. (4) is introduced briefly here. By assuming that a solution x(τ, ε) for Eq. (4) exists and is unique, we express it in the form as where x g is the general solution of the homogeneous equation and x p is a particular solution of the inhomogeneous equation. As in Ref. [33], the solution of the homogeneous equation can be obtained as where functions x 1 and x 2 denoted as form the fundamental system of solutions of the homogeneous differential equation. Besides, c 1 and c 2 are arbitrary constants. Further, a particular solution of the inhomogeneous equation can be obtained as In this work, different from Ref. [33]. where c 1 and c 2 are determined using the initial displacement and velocity, once the displacements x s and x e are given, the coefficients c 1 and c 2 can be determined as

Parameter identification for LTV dynamic systems
With the LTV dynamic system model given, as shown in Eq. (1), the objective of the present research is to identify the unknown stiffness k(t).
If the displacement, velocity, acceleration responses and external force can be measured directly and accurately, the identification process would be fairly simple and clear, which can be described by linear equations as in which t is the sampling time interval. However, there are two main problems with this identification strategy, which make it usually impossible to estimate the unknown parameter successfully.
(1) Eq. (12) is too sensitive to measurement noise. At each sampling time instant, even small measurement errors could lead to a large deviation from the true value. (2) It is usually unrealistic to measure the displacement, velocity, and acceleration signals simultaneously.
Generally, only one type of the response signals is measured in dynamic testing.

"Short time linearly varying" assumption for LTV dynamic systems
To address the aforementioned Problem (1), the main idea is to use fewer variables to represent the stiffness k(t), which can reduce the number of the unknowns in Eq. (12) and improve the robustness of identification to measurement noise. In a slowly varying system as mentioned in Sect. 2, the derivative of the parameter is a small variable. Thus, in a time interval with a certain width, the changing rate of the real parameter is relatively low. As a result, the TV parameter of such a slowly varying system can be fitted by a piecewise linear function with a relatively high precision. Therefore, the "short time linearly varying" assumption, which has a high accuracy and good capability of capturing the local variations of the model parameter [19], is adopted in this research. Under this assumption, a series of relatively short time windows are considered and then the TV stiffness k(t) within the ith time window can be written as in which t s i and t e i are the starting and ending time points of the ith time window, k i and s i are the assumed initial value and changing rate of k(t), respectively. With Eq. (13) substituted into model (1), the dynamic equation of the system within the ith time window can be rewritten as Thus, the total number of the unknown variables within the ith time window has been reduced to 2 by this assumption. A piecewise linear function will represent the estimated TV stiffness k(t) during the entire time period.

Using WKB method for parameter identification of LTV dynamic systems
In practical engineering, as mentioned in Problem (2), it is unrealistic to measure all of the displacement, velocity and acceleration signals simultaneously. For example, if acceleration signal is directly measured, numerical integration should be applied to the measurement to obtain the required velocity and displacement data. However, the measurement noise along with the cumulative error in the numerical integration would lead to large errors in the velocity and displacement data.
In the present research, the WKB method is used for parameter identification of LTV dynamic systems. Considering Eq. (4), with the "short time linearly varying" assumption adopted, the estimated k(τ ) within the ith short time window can be written as a linear function with respect to τ : where T s i and T e i are the starting and ending time points of the ith time window with respect to τ . Then, within the ith time window the considered dynamic equation can be written as Hence, in this short time window, the displacement response of Eq. (16) can be obtained by the WKB method.
For the sake of brevity, the derivation is given in Appendix 1 in detail. A good estimate of K i k i s i T within the ith time window can be obtained by finding the minimum of the Mean Absolute Error (MAE) [34][35][36], i.e.
where N i is the number of sampling points in the ith time window. It is worth noting that measured velocity or acceleration data can also be used in the identification process. In the ith time window, the velocity and acceleration responses of the considered dynamic system, v(τ ) and a(τ ), can be obtained using the analytical expression for the displacement response. The corresponding derivation is given in Appendix 1 in detail.Similar to x(τ ), v(τ ) and a(τ ) are also nonlinear functions with respect to k i and s i . To estimate K i with the measured velocityṽ(τ ) or acceleration signalsã(τ ) respectively, the following nonlinear optimization problems need to be solved: or Thus, only one type of response signal needs to be measured. No numerical integration or differentiation is applied to the measured signal. The WKB method can help radically avoid the errors caused by the truncation or cumulative error and the magnified noise in numerical calculation.

New LTV system identification based on WKB method
In this section, the application of the identification method combining the WKB method and "short time linearly varying" assumption will be realized for the LTV dynamic system.

Modification of WKB approximate solution in identification
As given in Appendix 1, in the forward problem for the LTV dynamic system, the coefficients c i,1 , c i,2 are determined by the measured displacement data at the starting and ending time points. However, in noise environment, the solutions obtained by Eq. (A5) would be highly sensitive to the measurement errors ofx(T s i ) andx(T e i ), which are inevitable in practice. The errors in measuredx(T s i ) andx(T e i ) are inevitable in practice.
In order to solve c i,1 and c i,2 more accurately in noise environment, the following linear least square equation is given The singular value decomposition can be adopted to determine the least square solution for c i,1 and c i,2 . Compared with the conventional solution, the modification made in Eq. (20) can effectively improve the robustness to the measurement noise. The values of c i,1 and c i,2 are not determined by the measured displacement at a certain time point. An accidentally large measurement error at one time instant will not strongly affect the solution accuracy. Equation (20) can provide more accurate c i,1 and c i,2 to the WKB approximation, which plays an important role in the nonlinear optimization.

Search strategy for nonlinear optimization problem in short time windows
In the ith time window, x p i ,x i,1 and x i,2 at each sampling point in Eqs. (A2)-(A3) cannot be expressed by any unified functions with respect to τ because the integration in Eqs. (A1)-(A4) cannot be resolved analytically. Therefore, the nonlinear optimization schemes widely used in nonlinear least square problems, such as Gauss-Newton and Levenberg-Marquardt, cannot solve Eq. (17). Thus, a new search strategy is proposed to solve it in this section. According to the derivation in Appendix 1, it can be obtained that Then, the objective is to find a k i and s i that can minimize the value of (21). Understandably, in the ith time window, the closer the line determined by k i and s i is to the true variation of k(τ ), the smaller is the value of (21), and vice versa. In this work, we define the percentage ratio of the Mean Absolute Error (MAE) with (k i , s i ) as R(k i , s i ) is used to visually measure the precision of the model output for the chosen k i and s i . For example, consider a LTV dynamic system: c 0.1, k 5 + 1.7 cos(0.5πεt), f (t) cos(10πεt), ε 0.1. Within the time window [1 s, 2 s], R(k i , s i ) with respect to k i and s i is shown in Fig. 1, which is a smooth surface. When the initial value of the supposed k(τ ) is fixed at 5, the dependency of R(5, s i ) on s i is shown in Fig. 2. Similarly, with the slope fixed at −0.2, the dependency of R(k i , −0.2) on k i is shown in Fig. 3. Both the curves given in Figs. 2 and 3 have only one extreme value.
In the search regions for k i and s i , [k min , k max ] and [s min , s max ] respectively, we use κ 1 , κ 2 , . . . , κ r and σ 1 , σ 2 , . . . , σ r that divide the search regions into a certain number of equal parts: First, R(κ 1 , σ 1 ), R(κ 1 , σ 2 ),. . ., R(κ 1 , σ r ) are calculated respectively and the smallest one R 1 (κ 1 ) R(κ 1 , σ q ) with respect to σ q is determined. According to Fig. 2 and the discussion above, the local optimal solution for s i when k i is fixed at κ 1 is in the interval σ q−1 , σ q+1 (or [s min , σ 2 ] if q 1, σ r −1 , s max if q r). Then, this interval is considered to be the new search region, and divided into equal parts. If the calculations and comparisons above are performed for the second round, the local optimal solution of s i at κ 1 can be determined in a smaller interval. By repeating this process j times until the latest calculated R j (κ 1 ) satisfies where the desired tolerance ρ is a small value, the smallest R for κ 1 can be determined as Similarly, R smallest can also be determined with the search process above. Among the calculated smallest Rs with respect to different initial values, the smallest one R smallest κ p can be determined. Then, the interval around κ p , κ p−1 , κ p+1 (or [k min , κ 2 ] if p 1, κ r −1 , k max if p r), is considered to be the new search region for k i . By repeating the search process introduced above, the optimal solution for k i and s i can be obtained.

An illustrative example: linear dynamic system with TV stiffness
A linear dynamic system with TV stiffness is considered to demonstrate the new proposed identification method with: c 0.1, k 5 + 1.7 cos(0.5πεt), ε 0.1, f (t) (1 + 0.05εt) cos(10πεt). The system is simulated by the fourth order Runge-Kutta method for 20 s with the response data sampled every 0.001 s. With the initial state x(0) x (0) 0, the simulated displacement, velocity and acceleration responses are regarded as the measured signal for the following identification process.
To verify the proposed method in this section, the traditional identification method proposed in Ref. [19] is selected as the baseline method. In the traditional method, under the "short time linearly varying" assumption, the TV stiffness is identified from the following linear least square equation where all the response signals x (τ ), x (τ ) and x(τ ) should be known in advance. The main idea of the traditional method is to solve the least square problem (26) using singular value decomposition.

Identification for numerical LTV system
In the identification process, the width of the time window is 1 s and the entire time period is divided as shown in Table 1. First, the simulated displacement response is taken as the measured signal. To evaluate the performance of the proposed method, the identification is carried out in the ideal condition without measurement noise, and also at different levels of white Gaussian noise (signal-to-noise ratios (SNRs) equal to 20, 15 and 10 dB), respectively.  The identification results for the entire time period in the four cases are shown in Fig. 4. In order to verify the method's advantages, the identification results obtained with the traditional method are also given for comparison. The values of MAE of k(τ ) in the non-dimensional form, defined as wherek(τ ) denotes the estimate of k(τ ) and N is the number of samples of the measured signal, are given in Table 2. According to Fig. 4 and Table 2, the identification results by the new method are generally closer to the real variation of k(τ ). In the case of no measurement noise, the two methods have similar identification precision. With the increased level of measurement noise, the new method still maintains relatively high identification accuracy, while the traditional method leads to much larger errors. The identification results are even out of the true variation range of the real stiffness, 3.3 to 6.7.
The main factor causing the errors in the traditional identification process is that the required velocity and acceleration are obtained with numerical differentiation. Thus, the truncation errors can be introduced  to the velocity and acceleration signals. Moreover, the measurement noise in the displacement signal will be magnified severely in the numerical differentiation. The effect of the measurement noise on the new method is limited. It exhibits good robustness to measurement noise.
The computation time of the proposed method and the traditional method for the considered case are 2.5 s and 0.3 s, respectively. The new method consumes more computation time, which is due to the optimization search algorithm. In the traditional method, the unknown parameters are obtained from a set of least square equations. The singular value decomposition can give the solution directly.
The identification results using the velocity and acceleration signals with different levels of measurement noise are shown in Figs. 5 and 6. The MAE of the identification results are given in Table 3. Compared with those using the displacement signal, the identification results have similar precision at different levels of measurement noise. The identification accuracy and robustness of the proposed method are thus also verified for velocity and acceleration measurement.

Effect of window size on identification accuracy
To investigate the effect of the window size on the identification results, two values, 0.5 s and 2 s, are also tested. The MAE of the identification results with different noise levels are shown in Table 4. For comparison, the MAE as in Table 2 for 1 s are also shown in Table 4. It is found that for the case of no measurement noise 0.5 s gives the best identification precision, MAE 0.56%. Increasing the window size will reduce the precision. For the cases of relatively low levels of measurement noise (SNR 20 dB and SNR 15 dB), 1 s gives the best identification precision, MAE 0.97% and 1.51%, respectively. While for the cases of high level of measurement noise (SNR 10 dB), 2 s gives the best identification precision. Decreasing the window size will reduce the precision. In the present numerical example, smaller time window has better accuracy when the noise level is low, and larger time window has better robustness when the noise level increases.
The errors with different window sizes are all acceptable. This is only a qualitative investigation based on a specific example. To draw a general conclusion about it, a rigorous theoretical proof is needed in further study.

Stability analysis of the proposed identification method
To verify the stability of the new identification method for systems with more rapidly changing parameters, k 5 + 1.7 cos(απεt), ε 0.1, α 0.5, 1, 2, 4 are considered respectively in this section. We set the window size to 0.5 s in this section. To focus on the stability with respect to the changing rate of k, the displacement signal with no measurement noise is taken. The MAE of the identification results for 0 to 20 s compared with those by the traditional method are shown in Table 5. Compared with the traditional method, the proposed method has much better accuracy when α increases. The main reason for this is that numerical differentiation or integration of the measurement signals is not involved in the identification process. Thus, the corresponding errors present in the traditional method can be avoided.

Conclusions
In this paper, a new LTV system identification method combining the WKB approximate solution and the "short time linearly varying" assumption is proposed, in which the TV stiffness is identified in short time windows with an optimization search algorithm. Numerical simulation is used to validate the proposed method. The identification results for cases of different levels of measurement noise are studied. The values of MAE are calculated to assess the estimation accuracy. The traditional method is also employed as a comparison to test the accuracy, robustness and stability of the new proposed method. The simulation results show that the developed search algorithm allows to estimate the TV stiffness with high precision, even when the measurement noise increases.
This work introduces the WKB approach to the LTV system identification for the first time. The new identification method has two distinct advantages as compared to the traditional one. First, only one type of signal (displacement, velocity or acceleration) needs to be measured in the identification process, which makes the practical application more convenient. Second, since the numerical differentiation and integration of the measured signal are not involved in the identification process, the truncation and cumulative errors can be avoided. This helps improve the robustness to the measurement noise and stability to the variation of unknown TV parameters.
It is worth noting that the modal decoupling of multiple degree of freedom (MDOF) system is not discussed in this paper. To experimentally validate the proposed method on practical structures, the modal decoupling of TV system should also be included, which is an important topic worth investigation in the future. Another limitation of the new proposed method lies in its computational efficiency. As discussed in Sect. 4.3, though the algorithm is well designed, it is somewhat time-consuming to find the optimal solution. However, this computational problem is less critical if a high-performance computing resources are available.