#### Abstract

This paper proposes a characteristic time expansion method (CTEM) for estimating nonlinear restoring forces. Because noisy data and numerical instability are the main causes of numerical developing problems in an inverse field, a polynomial to identify restoring forces is usually adopted to eliminate these problems. However, results of the way doing are undesirable for a high order of polynomial. To overcome this difficulty, the characteristic length (CL) is introduced into the power series, and a natural regularization technique is applied to ensure numerical stability and determine the existence of a solution. As compared to previous solutions presented in other researches, the proposed method is a desirable and accurate solver for the problem of restoring the force in the inverse vibration problems.

#### 1. Introduction

Identification of nonlinear dynamical system is a kind of inverse problems and is usually encountered in engineering applications. For instance, to specify the parameters of dynamical systems is necessary in optimal processes. It is important to analyze and determine the parameters of the system using experimental testing and numerical methods. However, uses of these methods might lead to some challenging problems in the structural mechanic field because a small measurement error can cause large errors in the results of parameters’ identification.

To overcome these inverse problems, some solvers in the literature were proposed by using numerical techniques and experimental testing. For instance, various publications in [1–10] recommended uses of the damping coefficient, stiffness, and external force for solving the inverse problems. Mode shape, frequency, displacement, and velocity at different times can also be used to estimate these properties successfully [11]. Huang [12] has employed a conjugate gradient method (CGM) to solve the nonlinear inverse vibration problem for the estimation of the time-dependent stiffness coefficient. Recently, Liu [13, 14] has developed a Lie-group shooting method to study the inverse vibration problem for estimating the time-dependent damping and stiffness coefficients and simultaneously derived a closed-form solution to estimate the parameters.

A complete review on the developments of useful methods for the realm of nonlinear system identification can be found in [14]. Reference [15] also proposed the idea of a force state mapping method which is a simple procedure and allows a direct identification of the restoring force for nonlinear mechanical systems. This idea was further extended in [16–18]. Recently, Namdeo and Manohar [19] have modified the force state mapping technique with two alternative schemes of functional representation: reproducing kernel particle method and kriging technique and estimating the parameters of nonlinear system from measured time histories of response under known excitations. Although this numerical method exhibits the capability to reproduce polynomials of specified order and has been applied to mechanical experiments, yet how to ensure numerical stability and avoid noise disturbance are not reported.

The aim of this paper is to develop a simple, multistep regularization algorithm with easy numerical implementation and versatility. A simple polynomial interpolation can be considered as a fit for the time history of displacement response under known excitations. Theoretically, if the order of the polynomial interpolation is large enough, the approximation displacement response can be as accurate as possible. However, in practical applications, it is not a good method to easily increase the order to obtain a highly accurate approximation. In fact, doing so will lead to a highly ill-posed matrix with a high-order function (the Vandermonde matrix), which has been described in [20]. To resolve the ill-posed linear problems for the Vandermonde system, Beckermann [21] and Li [22] claimed that the optimal condition number of the Vandermonde matrix could be expected. Nevertheless, the results of Beckermann [21] and Li [22] show that in the best possible cases, the condition numbers of the Vandermonde matrices still grow exponentially with the order of the interplant polynomial. Because of this, no one is interpolating the data by the high-order polynomials in the usual bases but rather in the Chebyshev polynomials. Hence, how to alleviate the ill-poseness owing to high-order function becomes one of the main tasks in this paper. First of all, we introduce the characteristic length (CL) of computational time into the high-order polynomial expansion to relieve the ill-conditioning of the resulting coefficient matrix of the polynomial expansion and then ensure numerical stability. This concept was first proposed to deal with the Laplace equation using a physical quantity [13, 14, 23, 24]. Recently, the CL has been successfully extended to deal with the Laplace equation and sloshing wave problems [25–27]. Although the CL can enhance the numerical accuracy for solving ill-posed linear matrix, it cannot avoid the effect of measured errors for parameters identification problems. Therefore, how to overcome the instability of the mathematical procedure is quite important. In addition, a small disturbance of measured data has to be considered in the numerical algorithm because they could cause an error identification of the parameter. In order to identify an accurate and stable solution for longer time scales, some special techniques have been used, including of the singular value decomposition (SVD), the SVD with a regularization parameter determined by the L-curve method, and sensitivity analysis. Despite these efforts, the stability problem remains unresolved. To thoroughly overcome these difficulties, this paper further adopts the CL combined with the natural regularization technique [28] to track ill-posed linear problems in numerical procedures. One advantage of this regularization method is that it can determine whether a solution exists for a linear system with the noisy level.

Apart from the current section, Section 2 describes the mathematical formulation of the characteristic time expansion method and introduces the numerical procedure of the matrix CGM. Section 3 gives several numerical examples, including of Duffing’s oscillator, Duffing’s oscillator with negative linear stiffness, Van der Pol’s oscillator, Bouc-Wen class model, and the seat model, to compare results of our method with the analytical solutions. Finally, some concrete conclusions are drawn in Section 4.

#### 2. Basic Formulation

A second-order ordinary differential equation (ODE) for the equation of motion is expressed as where represents the displacement of a system response; and are the external excitation and restoring force, respectively. In order to obtain , a trivial rearrangement of (1) gives Here, can be obtained if the quantities, and , on the right-hand side of (2) are known. In general, it is easier to measure the displacement at some discrete sampling times than to directly measure velocities and accelerations. Therefore, if is denoted as the measured displacement, the differentiation of displacements can be expressed as follows: This is however a set of index-three differential algebraic equations (DAEs) [29], which is hard to solve because the amplification of small errors and perturbations in the displacement cause severe numerically ill-posed condition.

##### 2.1. The Characteristic Time Expansion Method

The polynomial interpolation is defined as the interpolation of a given set of data by a polynomial. In other words, given some data points, the aim is to find a polynomial which exactly goes through these points of data.

According to (3), the displacement can be expressed as a polynomial expansion: where denotes each discrete time, denotes the displacement at each time, denotes the final time, and denotes the unknown coefficient. In many engineering applications, one wants to interpolate the data as accurate as possible. But this is limited by the interpolation of data with ()-order polynomials, where the resultant Vandermonde matrices are highly ill-conditioned as measured by the Lebesgue constant .

In this study, we introduce the characteristic length (CL) into (6) and express as follows where denotes the CL. Differentiation of (7) yields velocity and acceleration and they are expressed as follows: The polynomial expansion in (7) can be used to describe the displacement of a system. Hence, (7) can be expressed as a linear equation system with : We denote the above equation by where is the vector of unknown coefficients.

##### 2.2. The Matrix CGM for Ill-Posed Linear System

When a matrix is ill-posed and measured data contains noisy disturbances, it is difficult to ensure the stability of the system using the conventional regularization techniques. Therefore, Liu et al. [28] proposed a natural regularization method, which proves that a solution exists when ill-posed matrix and noisy disturbances occur. This method can be described by the following matrix equation: If is the inversion of , numerically, becomes a left-inversion of . Then we have Let Given , say , can be directly obtained because is a given matrix. Hence, we have When (11) and (14) are combined together, they create an over-determined system to calculate . The over-determined system can be written as where is an matrix with . Multiplying (14) by yields an matrix equation: Besides the primal system shown in (10), we need to solve the dual system with Applying the operators in (17) to and utilizing the above equation, that is, , we can obtain where .

Replacing the in (19) by , we have a similar equation for the primal system in (10) where .

Finally, when of (20) is calculated by the CGM, the restoring force, velocity and acceleration can be obtained from (2) and (8).

#### 3. Numerical Examples

*Example 1. *In this case, we consider a Duffing oscillator [29] and a second-order ODE to describe the forced vibration of a nonlinear structure by
where the parameters are fixed as , , and . The restoring force can be expressed as follows:
In order to identify the restoring force as a function of , a monotonic function of is required. In this example, is used to obtain the external force and is given by
To test the stability of the numerical method, the order of the polynomial and computational time are increased. The restoring force in the initial and final time changed very rapidly. To understand the CL effect, , , and are fixed. The maximum estimation error of with different CLs, shown in Figure 1, is smaller than .

It can be seen that including the CL into this case is efficient to overcome an ill-posed matrix. Furthermore, by fixing , the exact solutions for velocity and acceleration can be determined. The numerical results are shown in Figures 2, 3, and 4. According to the numerical results, the maximum estimation errors of are found to be smaller than . We can find that applying the CL and matrix regularization method with the CGM can provide highly stable and accurate solutions. In order to further test the stability of the present method under different noise levels, we also consider as an input into the estimation equations, where is a random number in , and is a noise level. With different noise levels %, 3%, 5%, and 10%, the computed profile of restoring forces is shown in Figure 5. Figure 5 also shows that the maximum estimated errors of are smaller than with noisy disturbances. We can find that the CL can effectively overcome numerical instability under the effect of the high-order function and large noise disturbances. Hence, we can see that the present method has a highly numerical accuracy and stability under the effect of the high-order function and large noise.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

*Example 2. *In this case, the Van der Pol oscillator is one of the nonlinear benchmark problem, and is given by
In this equation, is given by , and then, the external force can be obtained as
In this calculation, by fixing , , and , the numerical accuracy and stability of different parameters can be tested, including , and , , respectively. The maximum estimation error of shown in Figures 6 and 7 are smaller than.

**(a)**

**(b)**

**(a)**

**(b)**

To test the numerical stability of increasing the computational time by 10 seconds, the parameters are fixed as , and , , respectively. The maximum estimation errors of , shown in Figures 8 and 9, are smaller than . From the numerical solutions in Figures 6–9, it shows that the present method can keep the same numerical accuracy with the increase of the CL when the computational time increases.

**(a)**

**(b)**

**(a)**

**(b)**

This example demonstrates the results of fixing the parameters , , and under different noise levels with %, 3%, 5%, and 10%. The computed profile of is plotted in Figure 10. Figure 10(a) compares the restoring force with exact one, and the maximum estimation error of shown in Figure 10(b) is smaller than . From numerical result, we can find that the maximum error occurs because the displacement is equal to zero. That is, this present method can overcome the effect of the high-order function and large noise simultaneously. Therefore, it is found that the proposed method is accurate especially when noisy disturbances are encountered.

**(a)**

**(b)**

*Example 3. *The Bouc-Wen class model is one of the most widely used to efficiently describe smooth hysteretic behavior in engineering application. For a structural element described by the Bouc-Wen classical model, the restoring force is written as
where is a post- and preyield stiffness ratio, denotes natural frequency, and is an auxiliary variable that represents inelastic behavior. The evolution of is determined by an auxiliary ordinary differential equation, which can be written in the form of
where denotes the derivative of with respect to time and and are parameters that control the scale and sharpness of the hysteresis loops, respectively. Parameters, and , are used to control the shape of the hysteresis loop. In order to estimate the velocity and the restoring force of the Bouc-Wen model, the group preserving scheme (GPS), which was proposed by Liu [30], is adopted. The restoring force obtained by the GPS is referred to as the exact restoring force. We consider a system that has the parameter values as , , , , , , , , , and ; the initial condition of is given as (0.0, 0.0, 0.1). The exactly computed profile ofis plotted in Figure 11.

To obtain using the present method, the parameters , , , and are fixed. The computed profile of at 40 to 42 seconds is plotted in Figure 12(a), and the maximum estimated error of , shown in Figure 12(b), is smaller than . Further, under a noise of %, the computed profile of at 40 to 42 seconds is plotted in Figure 13(a). The maximum estimated error of , shown in Figure 13(b), is smaller than . We see from Figures 12 and 13 that the maximum estimated error of still keep in the order of under a noise of%. That is, we can use the present method to achieve a more accurate and stable solution under a large noisy level.

**(a)**

**(b)**

**(a)**

**(b)**

*Example 4. *As in Example 3, we consider the viscosity damping effect into the Bouc-Wen classical model and estimate the restoring force described as
where is the viscosity damping ratio. In this example, the parameter values are fixed as , , , , , , , , , , , and ; the initial condition of is given as . The exact computed profile of is plotted in Figure 14.

In this case, we use the same solver with the same parameters of Example 3. The computed profile of at 5 to 7 seconds is plotted in Figure 15(a). The maximum estimated error of , shown in Figure 15(b), is smaller than . Again, with a noise of %, the computed profile of at 5 to 7 seconds is plotted in Figure 16(a). The maximum estimated error of , shown in Figure 16(b), is smaller than . It can be seen in Figures 15 and 16 that the maximum errors are smaller than under a noise of %. Therefore, we conclude that for the smooth hysteretic behavior by the Bouc-Wen classical model, the accurate and stable solutions in Examples 3 and 4 are available when the proposed method is adopted.

**(a)**

**(b)**

**(a)**

**(b)**

*Example 5. *In order to test the numerical stability of the CTEM used for the restoring force problem of discontinuous type, the vehicle seat problem is considered. When the seated human body is exposed to vertical vibration, the single-degree of freedom model can be used to describe its seat-person mathematical behavior, as shown in Figure 17, and is given by [31]
The parameters of the discontinuous typed vehicle seat model are given as , , , , , and. Here the external force is given by , and the parameters are given by , , , and , respectively. The computed profile of by and 201 is shown in Figure 18(a). The maximum estimated error of , shown in Figures 18(b) and 18(c), is smaller than . Numerical results show that this present method does not exhibit the numerical oscillation (Gibb’s phenomenon) when the high-order function is used. Hence, this example demonstrates that the present method has a high level of accuracy and stability for the restoring force problem of discontinuous type.

**(a)**

**(b)**

**(c)**

#### 4. Conclusions

In nonlinear mechanical system analysis, the inverse vibration problem is difficult to solve under the measured data with noise. This paper has successfully combined the CTEM with a natural regularization algorithm to determine the unknown restoring force. Due to inclusion of the CL to retain high accuracy and stability, the approximation method can avoid the numerical instability caused by a high-order polynomial function. In addition, when the measured data is contaminated by a large noise, the errors can be controlled by utilizing a natural regularization technique and increasing the CL. In summary, the presented method is an effective and convenient approach to solve the inverse vibration problems.