An UKF‐based nonlinear system identification method using interpolation models and backward integration

In this paper, a novel identification method for nonlinear systems is proposed. This method utilizes linear interpolation models to describe the nonlinear forces of the physical models, and the unscented Kalman filter (UKF) method is adopted for the task of nonlinear identification. With the help of a linear interpolation algorithm, the proposed method requires little prior knowledge of the form of the nonlinear stiffness. Therefore, this method takes advantage of both the independence of the linear interpolation points and the inherent mathematical properties of the UKF. The UKF method is also modified to better fit the needs of parameter identification. To further emphasize parameter identification, backward integration and observations of the previous states are used. Two numerical simulations of the nonlinear elastic stiffness and Bouc–Wen hysteresis are conducted to show the flexibility and efficiency of this method. In these 2 examples, the observation signals are generated by analytic models, and the identifications are conducted with a linear interpolation model.

However, they are all built to better imitate physical behaviors, not for system identification. When these models are used in system identification, the identified systems are the projections of the real system based on the selected model basis.
Take the polynomial stiffness model as an example. Before parameter identification starts, the proper order of the polynomial should be determined. However, this identification is not an easy task, and many studies have focused on this topic. The least absolute shrinkage selection operator was introduced by Bonin, et al. [15] Korenberg, et al. developed a new effective method called the forward regression orthogonal estimator based on a forward regression procedure. [16] A method based on the Expectation Maximization algorithm was presented by Baldacchino, et al. [17] To some extent, these methods can help select the proper order of the polynomial; however, their computational cost is high.
The linear interpolation model has the advantage of intuitively describing the shape of the nonlinearity. It is purely linear between the interpolation points, and engineers have no need to select the model basis beforehand. These engineers need only estimate the number and position of the interpolation points. This work can be accomplished by gradually increasing the number of points, visually and iteratively. Moreover, this model is sufficiently accurate to approximate the nonlinear behavior in a civil engineering problem. The only challenging problem is the large number of parameters to be identified with the increasing number of interpolation point. The Bayesian class method shows a clever solution to this problem.
In Bayesian class methods, parameters are seen as random variables and are considered to be distributed in a range. With more parameters being allowed to change, less accuracy can be achieved for parameter identification. A good solution is that, in some steps, only a subset of all the parameters is updated until all the parameters are identified. This method reduces the number of related parameters and will solve this problem. In Section 4, a class of models based on linear interpolation is proposed. These models are piecewise linear, which largely reduces the nonlinearity of the mathematical expression. In this model, a maximum of three parameters are identified in each step. This model describes the shape of the nonlinear model, namely, nonlinear elasticity or hysteresis. Engineers only need to determine the nonlinear form of the target system before identification.
Considering Bayesian methods, the unscented Kalman filter (UKF) is a good solver for the nonlinear parameters identification problems, having the merits of lower computational costs, a sound theoretical foundation, and strong effect. Traditionally, researchers extend the use of the UKF to parameter identification by treating parameters as states of the system. [18,19] There are many kinds of modified UKF used for parameter identification. For example, in the iterative UKF, researchers intend to explore more information in a single step by iteratively correcting the state and parameters in one step. [20,21] The states of the system and parameters are estimated separately in the dual UKF method. [22,23] There are many adaptive UKF methods that adaptively estimate the covariance of the system noise and observation noise. [24,25] For some methods, the spread of the sigma points is adapted to cover the nonlinearity of the system. [26] S. S. Bisht successfully used the UKF method to track sudden changes of the stiffness in the multiple degree of freedom structure. [25] The UKF method was applied by M. S. Miah, et al. to identify the states and parameters of linear quadratic regulator control, and this control scheme is supported by experimental validation. The UKF identification method was also used for online multiple degree of freedom nonlinear structure identification by Thaleia Kontoroupi and A. W. Smyth. [27] R. Astroza, et al. identified the parameters of the finite element model using the UKF method and the Opensees software. [3] In most of these methods, the observation computed with the states of the current step (CO, short for the current observation) is used to update the states. In this paper, the modified UKF methods of these class are called standard UKF methods. A new version of UKF modified observations are developed and are formed by the CO and COs obtained from the former steps (CPO, short for current and previous observation). These are called backward UKF methods.
For the standard UKF methods, because some parameters are not directly presented in the calculation of the observations, the performance of tracking the observations outperforms that of the identification of the parameters. A. W. Smyth, et al. found that, in some cases, misleading result may be obtained because it is common that accurate tracking of measurements is achieved without complete parameter convergence. [28] This method uses backward integration to magnify the difference between the currently identified system and the target system, so that the goal of identifying the parameters consistently agrees with the tracking of observations. With this modification, the identification of the parameters is further emphasized, and a better identification is acquired.
In this paper, structures with single degrees of freedom (SDOF) are considered. The linear interpolation method is used to describe the nonlinear part of the nonlinear system, and backward integration is used to improve the performance of the UKF in identifying the parameters of the interpolation model. The structure of this paper is as follows. In Section 2, the algorithm of the scaled UKF as an example of the standard UKF methods is briefly reviewed. The insufficiency of the standard UKF for parameter identification is stated. In Section 3, the introduced method with the interpolation model, linkage strategy, and backward integration are developed. In Section 4, two numerical SDOF nonlinear structures with nonlinear elastic stiffness and hysteresis, respectively, are studied. The observation signals of these SDOF structures are generated by analytic models, whereas the identification is carried out with the linear interpolation model, which demonstrate the effectiveness and wide applicability of the proposed method.

| Unscented transform
The UKF relies on the unscented transform to propagate the mean and covariance of a random variable x through a nonlinear transform y = f(x). In the unscented transform, sigma points are used to estimate the mean and covariance of y given the mean and covariance of x. Sigma points are usually generated by the mean and covariance of x, as shown in Equation 2. The scaled UKF is an improved version of the original UKF, belongs to the standard UKF class and was first studied by Simon J. Julier. [29][30][31] When x has L elements, Julier constructed 2L + 1 sigma points to estimate the transformed mean y and covariance P y . The sigma points are generated as follows: (1) where α ∈ (0, 1] is a scaling parameter, which was set as 0.01 in this article. ffi p calculates the square root of matrix, which is defined as P ¼ ffiffiffi P p T ffiffiffi P p . [32] { } i , returning the i-th column of a matrix. After 2L points are obtained, the mean and covariance after a nonlinear transform can be estimated as where β was set as 2, and the distribution of x is assumed to be Gaussian. The superscripts (m) and (c) represent the weighting coefficients when estimating the mean and covariance. The superscript T denotes the transposition of a matrix. For non-Gaussian variables, the accuracy of third and higher order moments is determined by the choice of α and β. For details of the derivation and proof, please refer to Simon J. Julier's book. [29] 2.2 | UKF for a SDOF nonlinear stiffness system The procedure of the scaled UKF is presented in this subsection. Consider a discrete state space system as follows.
where x k is the state value in the k-th step with dimension L, y k is the observation of the system, w k is the vector of the additive system noise, v k is the additive observation noise, and F(·) is a nonlinear transform, and C is the observation matrix.
The procedure of the scaled UKF is summarized in Table 1.

| Linear interpolation model
In this section, a model based on the linear interpolation algorithm is introduced in detail. The UKF method utilizes cross covariances between the states and observations to approximate the differences between the current identified system and the true system. In the interpolation method, the points that are not adjacent to each other are independent; when one observation of the CPO is used, in each step, only two parameters are related. Consequently, in each updating step, only some of the parameters are updated. Moreover, the interpolation model is a piecewise linear model, which does not have severe nonlinearities. When applied for identification, this process greatly simplifies the complexity of a numerical model of a target system. The interpolation model describes only the shape of the stiffness and reflects the behavior of the nonlinear model intuitively. For other models, the estimation results are usually the projection of the real model on the space of the selected basis space. The selection of the basis, taking the order of the polynomial model as an example, is vital, and this range is usually difficult to choose. A broad range results in tedious and complex expressions and severe nonlinearity, although a small one has the possibility of losing information. The linear interpolation model is used for moderately varying nonlinear behaviors. For nonlinear behaviors that change severely, the interpolation model and other models suffer identical problems. This shows the advantage of the linear interpolation model over the traditional analytical model.

| Interpolation model for nonlinear elastic behavior
increasing the number of points. This model can easily imitate the asymmetric behavior that is very common in real applications and that is difficult for traditional models to address, whereas the traditional models exhibit tedious and complex forms. For symmetric models, the number of the model parameters can be reduced by half.
The mathematical expression of this model is where g(x i ) is the ordinate of the interpolation point at x i . x i is the abscissa of the interpolation point, x is the input, and g(x) is the output of the model. In this paper, CPO is used. In Figure 2, the star stands for the CO of the system, and the asterisks stand for the former time history of the system observations; the other markers have the same meanings as in Figure 1. The length of the CPO is selected such that the corresponding related interval in the interpolation model is within two interpolation points. A longer length increases the effect of backward integration but also increases the computational effort, and vice versa.
These two points are marked by enlarged circles, and the length of the CPO is chosen as three in Figure 2. There are, at most, three parameters related to such a choice of the CPO given the condition that the head of the CPO has just traveled to the next interval, or that the tail of the CPO has not entirely left the latest interval. Moreover, in the UKF method, a cross covariance matrix of the states is used to calculate the direction of the updating step. In the cross-covariance matrix, the value corresponding to the states unrelated to the observations are zero. Accordingly, the Kalman gain of the parameters of the unrelated points are zero, which means that they are not updated in the current step. As a result, the states of the filter can be separated into four groups: the states of the system, the parameters that do not belong to the interpolation model, the parameters that are related to the CPO in the interpolation model, and the parameters that are not related to the CPO in the interpolation model. It is obvious that only the first three groups are actually useful for identification. There are many kinds of mathematical model for hysteretic behavior, such as the Iwan model, [33] Prandtl-Ishlinskii model, [34] Preisach model, [35] Dahl model, [36] Duhem model, [37] and Bouc-Wen model. [12] Each model has its advantages and can imitate different kinds of hysteretic behaviors. The Bouc-Wen model has been applied in many fields. It has a simple expression that is easy to use and can describe a large amount of hysteretic phenomena. The analytical expression of the Bouc-Wen model is In this model, A, β, γ, and n are the parameters of this model; z is the restoring force, which is the output, and _ x stands for the velocity, which is the input of this model. n is an important parameter that increases the complexity of basis function of the Bouc-Wen model, making n most difficult to decide. The interpolation model can solve this problem because it does not need such a complex expression. The interpolation hysteresis model is expressed as follows.

| Linkage strategy for zero initial values
For the identification methods based on Kalman filter class, a good initial value is very important. It is very difficult to set the initial values for the parameters of linear interpolation model. Experience from many numerical experiments shows that a zero initial value acquires good results when identifying the parameters of the elastic stiffness. However, when a hysteresis model is studied, zero initial values cannot work. Thus, a numerical trick called a linkage strategy is introduced in this subsection.
Take the linkage strategy of a hysteresis model as an example, as is described in Equations 11 to 13. Assume the numbers of interpolation points for both l 1 (s) and l 2 (s) are N i = N d = 7. The interpolation points can be divided into two groups, s 1 n i and s 2 n d , with n i ∈ [1, N i ] and n d ∈ [1, N d ]. l 1 (s) and s 1 n i are used when z _ x≥0, and l 2 (s) and s 2 n d are used when z _ x< 0. An index is used for each interpolation point to mark whether the input |z| has arrived at the interval formed by the corresponding interpolation point. For example, assume an interval is formed by s 1 n 1 ¼ 0 and s 1 n 2 ¼ 1; if the input |z | ∈ [0, 1] and z _ x≥0, then we state |z| has arrived at the interval formed by s 1 n 1 and s 1 n 2 . Initially, the indexes of the interpolation points are all set to be 1, which means that the input never arrived at any of the interval formed by these points. In addition, the value of the ordinates of the interpolation points is acquired by the rules described as follows.
When the input arrives at one interval, the indexes of the related two points are set to 0. The ordinates of these points will then be acquired by the UKF, and the ordinates of the other points whose indexes are 1 follow those of the interpolation points with index 0. The initial condition of the linkage strategy is shown in Figure 3.
• Condition 1: when excitation starts, the hysteretic force z starts to depart from zero. Assume that, at the beginning, z _ x ≥ 0. Then, as is shown in the right side of Figure 4, |z| has arrived at the interval s 1 1 ; s 1 2 Â Ã and is increasing; then, the indexes of the interpolation points s 1 1 and s 1 2 are set to be 0. This means that the ordinate of points s 1 1 and s 1 2 should be set as the value identified by the UKF, and the ordinate of the other interpolation points in the group s 1 n i will follow the ordinate of s 1 2 .
• Condition 2: as shown in Figure 5, the corresponding state has increased into the interval s     Figure 6 will be set to be the same as that of s 2 4 , and the ordinates of points to the right of s 2 5 will be set to be the same as that of s 2 5 .
The linkage strategy eliminates the need for the initial value estimation and increases the speed of identification. For the elastic interpolation model, according to the experience of the authors, the linkage strategy is also helpful but not necessary.

| Backward integration
This research focuses on the parameter identification of nonlinear dynamic systems with SDOF, the basic form of the target model is The structure is normalized to have a unit mass. In Equation 14, c is the equivalent viscous damping coefficient. f (x, history of x) represents the nonlinear force that is only related to current and the previous displacement. This force can either be elastic or hysteretic.
The CPO heightens the difference caused by the different parameters, which usually hides behind the system states. However, too many points increase the burden of calculation when computing the observation correlation matrix and the Kalman gain matrix and sometimes bring numerical errors into the estimation when the condition number of the observation correlation matrix is large. Therefore, a fixed length CPO is used, and the rule to choose this length are mentioned in the Subsection 3.1.1.
In the procedure of a typical Kalman filter method, the observations are estimated with the updated states, which contain the states of the system and the parameters to be identified. Some authors used the updated parameters to calculate the system time history from the zero initial conditions to the current step. This method has the hidden assumption that the parameters are constant, which is not always the case. The parameters could vary slowly with time. Accordingly, to the updated parameters with the newest observation and implement backward UKF method, the observations of the previous steps should be estimated with the current states and parameters. The backward integration is the time reverse process of the forward integration and will solve this problem. The derivation of the backward integration is as follows.
Assuming p(t c − t) = x(t), where t c stands for the starting time of the backward integration, and define m = t c − t. Then, the backward velocity is given as follows.
For the second derivative of x with respect to t, € x ¼ d dx dt =dt, we can derive the backward acceleration.
Substituting Equations 15 and 16 into Equation 14, the second order backward differential equation can be expressed as follows.
The expression of f (p, 'future' of p) needs further derivations when dealing with different nonlinear models. An example of the backward integration of a system with a nonlinear elastic stiffness will demonstrate the accuracy of backward integration.
The example is a nonlinear system: In this example, f (x) = k 1 x + k 3 x 3 , where k 1 and k 3 are the coefficients of the cubic nonlinear stiffness. The excitation signal € x g was a 1 Hz sine wave. The integration time step Δt was 0.001 s. The forward and backward integration was carried out for 1,000 steps, with the numerical integration method Runger-Kutta 4 (RK4). The last states of the forward integration, x(t c ) and − _ x t c ð Þ, are used as the initial conditions of the backward integration. In Figure 7, the results of the forward integration and backward integration are shown to coincide with each other. The error analysis in Figure 8 that subtracts the result of the backward integration from the forward integration shows that the error ratio is on the 10 −8 order, which can be negligible.
In the last example, the linear part of the system possesses a damping ratio of ξ ¼ c 2 ffiffiffiffi ffi k 1 p ¼ 0:03. The backward integration utilizes the states of the current time to estimate the states of the former time, which is different from forward integration. Obviously, when large damping is involved, the accuracy will decrease. However, the real system is accurate with no numerical error, and if we decrease the time step size, the accuracy of the backward integration will increase at the same time. In the example shown in Figure 9, the damping ratio is increased to 1, the time step size is decreased to 0.0001 and the results for 10,000 steps of forward and backward integration are compared. Figure 10 also shows the high accuracy of the case when the damping ratio is very large.

| Summary of the proposed method
In this section, the above methods will be integrated into the UKF algorithm. The SDOF nonlinear elastic stiffness system is selected as an example.
The target system is the same as in Section 2.2. g(x) represents the interpolation model of the nonlinear stiffness. The integrated model is where g(x) is the same as in Equation 9 .  The first order expression of Equation 19, considering the system noise, is The corresponding discrete system is with where θ = [c, x 1 , x 2 , …, x n ] T is the parameter vector that consists of the damping coefficient and the parameters of the interpolation model. ΒI:R n →R n l (or R k when k < n l ) is a backward integration operator mapping the current states and parameters to the former n l step time history of the observed variable. k is defined in Table 1, and n l is the length of the CPO. The symbol R in this mapping should be distinguished from the covariance matrix of the observation noise R. The other symbols have the same meanings as introduced in Section 2. Then, the procedure of the proposed method can be summarized in Table 2. In this table, only Steps 1, 2, 4, and 7 are different from the UKF algorithm introduced in Table 2.

| Model identification of a nonlinear elastic model
In this section, an SDOF nonlinear stiffness system is used as the target system, and the damping coefficient and the parameters of the interpolation model are identified. The target system is The initial conditions were x = 0, _ x ¼ 0, and n l , and the length of the CPO was set as 60. The sampling rate and the integration time step were 0.001 s; the RK4 integration method was used for both the forward integration and the backward integration. The corresponding state space expression of the nonlinear interpolation model is where The CPO had 60 points, which consist of the newest displacement and the recent 59 displacement history, and an additive noise level of 5% of the root mean square. The excitation signal is a 1 Hz sine wave.
The identified interpolation stiffness model is presented in Figure 11. In Figure 11, the identified model agrees very well with the true model. The interpolation points whose intervals have not been reached by state x stay in arbitrary places, which might be caused by false correlation values involved in numerical error when calculating the inverse of the observation covariance matrix. If the behavior outside the range of the states is considered, extrapolation could be used. FIGURE 11 Comparison of the identified stiffness with the true stiffness Figure 12 depicts the convergence process of the interpolation points. When the first appears in the interval formed by the two interpolation points, they move quickly to the right position and stay there for the rest of the identification process, which also verifies the presumption that the parameters that are not adjacent to each other are irrelevant.
As is shown in Figure 13, the damping coefficient changes as the parameters of interpolation model change and stays steady when the parameters of the interpolation model converge. The last identified value was not adopted as a result of    Figure 14 statistically demonstrates the steady feature of this variable. The multiplicity of the histogram is 1.13, which shows a high accuracy compared with the true parameter. The interval of the histogram is 0.01; therefore, the error is 0, which means that the error ratio is less than 1%.

| Model identification of a hysteresis model
In this section, a SDOF hysteretic system was chosen as the target system The initial conditions were x = 0, _ x ¼ 0, and z = 0. The length of backward integration was set as 60. The sampling rate and the integration time step were 0.001 s; the RK4 integration method was used for both the forward integration and the backward integration. The corresponding state space expression of the nonlinear interpolation model is Where s 1 is the ordinate of the interpolation points that contains 21 points with their horizontal coordinates ζ 1 evenly spaced over [0, 40], and s 2 has 8 points that are horizontally evenly spaced in [0, 40] with the abscissa vector ζ 2 . This choice was made because the increase of the hysteretic force was much slower than its decrease. Although the identification of the decreasing process is relatively rough, as demonstrated in Figure 17, the estimation is precise enough to represent the behavior of the hysteretic process. Because the hysteretic force is related with time, the backward time integration differential equation needs to be derived.
Substitute Equations 38 and 15 into Equations 31 and 34. The backward integration differential equations are written as follows. (39) The CPO was set as x and its 60 recent time history points, with the additive noise level of 5% root mean square. The excitation signal € x g was the combination of three sine waves with their frequencies at 1, 2, and 3 Hz, respectively. This choice was intentionally set to fully stimulate the characteristics of the hysteretic behavior. Figure 15 shows that, for the first 500 steps of the states tracking, the proposed method traces the states of the target system very well and shows a very quick convergence rate. A good tracking of the hysteretic force is the premise of a correct estimate of the hysteretic model. Figure 16 illustrates the convergence process of the parameters of the interpolation model. All the parameters of s 1 leave zero at the same time and quickly converge to the correct position. When the input of the interpolation model decreases, all the parameters of s 2 leave zero at the same time. This indicates that the linkage strategy works very well. With the vibration of the states, the parameters slowly converges to the correct position. The convergence rate is not as fast as that in Section 4.2. The reason for this is that the state z is itself a nonlinear first order system. The sensitivities of the parameters in the interpolation model with respect to the CPO are reduced by the uncertainty FIGURE 15 First 500 steps tracking the process of the displacement, velocity, hysteretic force of the system of the state z. Therefore, the parameters change much slower when moving towards the most feasible position than occurs in Section 4.2.
In Figure 17, the solid line represents the true input and output relationship. The vertical line stands for the jump between the two processes described in Equation 34. The circles describe the identification results of s 1 , and the asterisks  stand for the identification result of s 2 . The picture shows high coincidence between the estimated model and the true model. There are intervals that are never reached by the input. Their ordinates are assigned the same values as those of the nearest identified parameter. This phenomenon also indicates the strong performance of the linkage strategy. Figure 18 gives the estimation of the damping coefficient as the value 1.17. The true value is 1.13, and the error ratio is computed as 3.5%. This shows the accuracy and efficiency of this new identification method.

| CONCLUSION
In this article, an identification method based on an UKF method is proposed. An interpolation model is introduced in Section 3 to imitate the behaviors of nonlinear elastic stiffness and hysteresis. To integrate this model with the UKF identification method, a linkage strategy and backward integration are derived. This method involves using historical observation points to tell the difference between the true system and the currently identified system. In Section 4, another two numerical examples, which are a nonlinear elastic stiffness SDOF system and a hysteretic SDOF system, are adopted to verify the effectiveness of the proposed method. The parameters and states converge very fast with high accuracies.
This method has the following advantage as.
1. A complex function form for the physical behavior is not needed. The interpolation model used in this method is the linear interpolation functions that depict the shape of the physical behavior. 2. The proposed method does not need to choose the proper order of the model or to use a large number of computations to find this order. 3. The parameters of the interpolation model are not always related. Taking advantage of the inherent mathematical properties of the UKF method, only a small number of parameters are identified.