# Application of the Characteristic Time Expansion Method for Estimating Nonlinear Restoring Forces

Authors: Yung-Wei Chen
Journal: Journal of Applied Mathematics
 (2013)

---

## 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.


---

## Body

## 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: (1) reproducing kernel particle method and kriging technique and (2) 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, Section2 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(1)x¨+H(x,x˙)=P(t),
where x represents the displacement of a system response; P(t) and H(x,x˙) are the external excitation and restoring force, respectively. In order to obtain H, a trivial rearrangement of (1) gives
(2)H(x,x˙)=P(t)-x¨.
Here, H can be obtained if the quantities, P(t) and x¨, 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 x1(t)=g(t) is denoted as the measured displacement, the differentiation of displacements can be expressed as follows:
(3)x1(t)=g(t),(4)x˙1(t)=x2(t),(5)x˙2(t)=x3(t).
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:
(6)x(tk)=∑k=1nak(tk)n-1,0≤tk≤tf,
where tk denotes each discrete time, x(tk) denotes the displacement at each time, tf  denotes the final time, and ak 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 n data with (n-1)-order polynomials, where the resultant Vandermonde matrices are highly ill-conditioned as measured by the Lebesgue constant 2n/[e(n-1)lnn].In this study, we introduce the characteristic length (CL) into (6) and express as follows
(7)x(tk)=a0+∑k=1mak(tkT0)k,0≤tk≤tf,tf<T0,
where T0 denotes the CL. Differentiation of (7) yields velocity and acceleration and they are expressed as follows:
(8)x˙(tk)=∑k=1mkT0ak(tkT0)k-1,x¨(tk)=∑k=1mk(k-1)T02ak(tkT0)k-2.
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 n=m+1:
(9)[1t1T0(t1T0)2⋯(t1T0)m1t2T0(t2T0)2⋯(t2T0)m1t3T0(t3T0)2⋯(t3T0)m⋮⋮⋮⋱⋮1tmT0(tmT0)2⋯(tmT0)m][a0a1a2⋮am]=[x(t1)x(t2)x(t3)⋮x(tm)].
We denote the above equation by
(10)Rc=b1,
where c=[a0,a1,a2,…,am]T  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:
(11)RTUT=Im,thatis,(UR)T=Im.
If U is the inversion of R, numerically, U becomes a left-inversion of R. Then we have
(12)(RRT)UT=R.
Let
(13)RX0=y0.
Given X0, say X0=Ι=[1,…,1]T, y0 can be directly obtained because R is a given matrix. Hence, we have
(14)y0TUT=X0T,thatis,X0=Uy0.
When (11) and (14) are combined together, they create an over-determined system to calculate UT. The over-determined system can be written as
(15)BUT=[ImX0T],
where
(16)B:=[RTy0T]
is an n×m matrix with n=m+1. Multiplying (14) by BT yields an m×m matrix equation:
(17)[RRT+y0y0T]UT=R+y0X0T.
Besides the primal system shown in (10), we need to solve the dual system with
(18)RTy=b1.
Applying the operators in (17) to b1 and utilizing the above equation, that is, y=RTb1, we can obtain
(19)[RRT+y0y0T]y=Rb1+(X0·b1)y0,
where y0=RX0.Replacing theR in (19) by RT, we have a similar equation for the primal system in (10)
(20)[RTR+y0y0T]c=RTb1+(X0·b1)y0,
where y0=RTX0.Finally, whenc of (20) is calculated by the CGM, the restoring force, velocity and acceleration can be obtained from (2) and (8).


## 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:
(6)x(tk)=∑k=1nak(tk)n-1,0≤tk≤tf,
where tk denotes each discrete time, x(tk) denotes the displacement at each time, tf  denotes the final time, and ak 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 n data with (n-1)-order polynomials, where the resultant Vandermonde matrices are highly ill-conditioned as measured by the Lebesgue constant 2n/[e(n-1)lnn].In this study, we introduce the characteristic length (CL) into (6) and express as follows
(7)x(tk)=a0+∑k=1mak(tkT0)k,0≤tk≤tf,tf<T0,
where T0 denotes the CL. Differentiation of (7) yields velocity and acceleration and they are expressed as follows:
(8)x˙(tk)=∑k=1mkT0ak(tkT0)k-1,x¨(tk)=∑k=1mk(k-1)T02ak(tkT0)k-2.
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 n=m+1:
(9)[1t1T0(t1T0)2⋯(t1T0)m1t2T0(t2T0)2⋯(t2T0)m1t3T0(t3T0)2⋯(t3T0)m⋮⋮⋮⋱⋮1tmT0(tmT0)2⋯(tmT0)m][a0a1a2⋮am]=[x(t1)x(t2)x(t3)⋮x(tm)].
We denote the above equation by
(10)Rc=b1,
where c=[a0,a1,a2,…,am]T  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:
(11)RTUT=Im,thatis,(UR)T=Im.
If U is the inversion of R, numerically, U becomes a left-inversion of R. Then we have
(12)(RRT)UT=R.
Let
(13)RX0=y0.
Given X0, say X0=Ι=[1,…,1]T, y0 can be directly obtained because R is a given matrix. Hence, we have
(14)y0TUT=X0T,thatis,X0=Uy0.
When (11) and (14) are combined together, they create an over-determined system to calculate UT. The over-determined system can be written as
(15)BUT=[ImX0T],
where
(16)B:=[RTy0T]
is an n×m matrix with n=m+1. Multiplying (14) by BT yields an m×m matrix equation:
(17)[RRT+y0y0T]UT=R+y0X0T.
Besides the primal system shown in (10), we need to solve the dual system with
(18)RTy=b1.
Applying the operators in (17) to b1 and utilizing the above equation, that is, y=RTb1, we can obtain
(19)[RRT+y0y0T]y=Rb1+(X0·b1)y0,
where y0=RX0.Replacing theR in (19) by RT, we have a similar equation for the primal system in (10)
(20)[RTR+y0y0T]c=RTb1+(X0·b1)y0,
where y0=RTX0.Finally, whenc 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
(21)x¨+γx˙+βx+αx3=P(t),
where the parameters are fixed as α=1, β=-1, and γ=0.3. The restoring force can be expressed as follows:
(22)H(x)=x3-x.
In order to identify the restoring force H as a function of x, a monotonic function of t is required. In this example, x(t)=t2-8 is used to obtain the external force and is given by
(23)P(t)=(t2-8)3-t2+0.6t+10.
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, m=201, X0=I, and ε=1×10-16 are fixed. The maximum estimation error of H with different CLs, shown in Figure 1, is smaller than 10-6.Figure 1
For Example1, showing the error of estimation with different CLs.It can be seen that including the CL into this case is efficient to overcome an ill-posed matrix. Furthermore, by fixingT0=1200, 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 H are found to be smaller than 6×10-10. 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
(24)x^i=xi+σR(i)
as an input into the estimation equations, where R(i) is a random number in [-1,1], and σ is a noise level. With different noise levels σ=1%, 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 H are smaller than 10-1 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.For Example1, (a) comparing the exact solution and numerical one for velocity and (b) displaying the estimated error of velocity.
(a)
(b)For Example1, (a) comparing the exact solution and numerical one for acceleration and (b) displaying the estimated error of acceleration.
(a)
(b)For Example1, (a) comparing the exact solution and numerical one for restoring force and (b) displaying the estimated error of restoring force.
(a)
(b)For Example1, (a) comparing estimated and exact restoring forces under different noise level and (b) displaying the error of estimation with σ=1%, 3%, 5%, and 10%, respectively.
(a)
(b)Example 2.
In this case, the Van der Pol oscillator is one of the nonlinear benchmark problem, andH(x,x˙) is given by
(25)H(x,x˙)=x+(x2-1)x˙.
In this equation, x is given by x(t)=t3/3-8, and then, the external force can be obtained as
(26)P(t)=(t33-8t)+[(t33-8t)2-1](t2-8)+2t.
In this calculation, by fixing ε=1×10-16, X0=I, and Tf=4, the numerical accuracy and stability of different parameters can be tested, including T0=101, m=100 and T0=2×105, m=500, respectively. The maximum estimation error of H shown in Figures 6 and 7 are smaller than10-8.For Example2, (a) comparing exact solution and numerical one for restoring force withT0=101 and m=100 under 4 seconds and (b) displaying the estimated error of restoring force.
(a)
(b)For Example2, (a) comparing exact solution and numerical one for restoring force with T0=2×105 and m=500 under 4 seconds and (b) displaying the estimated error of restoring force.
(a)
(b)To test the numerical stability of increasing the computational time by 10 seconds, the parameters are fixed asT0=120, m=100 and T0=4×105, m=500, respectively. The maximum estimation errors of H, shown in Figures 8 and 9, are smaller than 10-8. 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.For Example2, (a) comparing exact solution and numerical one for restoring force with T0=120 and m=100 under 10 seconds and (b) displaying the estimated error of restoring force.
(a)
(b)For Example2, (a) comparing exact solution and numerical one for restoring force with T0=4×105 and m=500 under 10 seconds and (b) displaying the estimated error of restoring force.
(a)
(b)This example demonstrates the results of fixing the parametersTf=10, T0=1.2×104, and m=201 under different noise levels with σ=1%, 3%, 5%, and 10%. The computed profile of H is plotted in Figure 10. Figure 10(a) compares the restoring force with exact one, and the maximum estimation error of H shown in Figure 10(b) is smaller than 10-1. 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.For Example2, (a) comparing estimated and exact restoring forces and (b) displaying the error of estimation with σ=1%, 3%, 5%, and 10%, respectively.
(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(27)H(x,x˙)=αwn2x+(1-α)wn2z,
where α is a post- and preyield stiffness ratio, wn denotes natural frequency, and z is an auxiliary variable that represents inelastic behavior. The evolution of z is determined by an auxiliary ordinary differential equation, which can be written in the form of
(28)z˙=Ax˙-βx˙|z|n-γ|x˙||z|n-1z,
where z˙ denotes the derivative of z with respect to time and A and n 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 A=0.5, β=-5.0, γ=5.0, n=1.4, α=0.4, wn=3.0, t0=0.0, tf=50.0, Δt=0.01, and P(t)=sin(t); the initial condition of (x,x˙,z˙) is given as (0.0, 0.0, 0.1). The exactly computed profile ofHis plotted in Figure 11.Figure 11
For Example3, restoring force surface based on GPS.To obtainH using the present method, the parameters T0=8, X0=I, m=201, and ε=1×10-14 are fixed. The computed profile of H at 40 to 42 seconds is plotted in Figure 12(a), and the maximum estimated error of H, shown in Figure 12(b), is smaller than 1.4×10-2. Further, under a noise of σ=5%, the computed profile of H  at 40 to 42 seconds is plotted in Figure 13(a). The maximum estimated error of H,  shown in Figure 13(b), is smaller than 2.5×10-2. We see from Figures 12 and 13 that the maximum estimated error of H still keep in the order of 10-2 under a noise ofσ=5%. That is, we can use the present method to achieve a more accurate and stable solution under a large noisy level.For Example3, (a) comparing estimated and exact restoring forces and (b) displaying the estimation error of restoring force.
(a)
(b)For Example3, (a) comparing estimated and exact restoring forces with σ=5% and (b) displaying the estimation error of restoring force.
(a)
(b)Example 4.
As in Example3, we consider the viscosity damping effect into the Bouc-Wen classical model and estimate the restoring force described as
(29)H(x,x˙)=2ξwnx˙+αwn2x+(1-α)wn2z,
where ξ is the viscosity damping ratio. In this example, the parameter values are fixed as A=0.8, β=4.0, γ=2.1, n=1.4, α=0.4, wn=3.0, ξ=0.15, t0=0.0, tf=10, Δt=0.01, X0=I, and P(t)=0.4sin(t); the initial condition of (x,x˙,z˙) is given as (0,0,0.1). The exact computed profile of H is plotted in Figure 14.Figure 14
For Example4, restoring force surface based on the GPS.In this case, we use the same solver with the same parameters of Example3. The computed profile of H at 5 to 7 seconds is plotted in Figure 15(a). The maximum estimated error of H, shown in Figure 15(b), is smaller than 2.5×10-3. Again, with a noise of σ=5%, the computed profile of H at 5 to 7 seconds is plotted in Figure 16(a). The maximum estimated error of H, shown in Figure 16(b), is smaller than 1×10-2. It can be seen in Figures 15 and 16 that the maximum errors are smaller than 10-2 under a noise of σ=5%. 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.For Example4, (a) comparing estimated and exact restoring forces and (b) displaying the estimation error of restoring force.
(a)
(b)For Example4, (a) comparing estimated and exact restoring forces with σ=5% and (b) displaying the estimation error of restoring force.
(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 Figure17, and is given by [31]
(30)Mx¨+c1x˙+c2|x˙|x˙+k11+k2|x|x=P(t).
The parameters of the discontinuous typed vehicle seat model are given as k1=48000Nm-1, k2=24000Nm-1, c1=300Nsm-1, c2=1500Nsm-1,  M=8kg,  andM1=42kg. Here the external force is given by P(t)=0.04cos(t), and the parameters are given by T0=11, tf=10, X0=0.001, and ε=1×10-14, respectively. The computed profile of H by m=51 and 201 is shown in Figure 18(a). The maximum estimated error of H, shown in Figures 18(b) and 18(c), is smaller than 5×10-3. 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.Figure 17
For Example5, seat-person model of single-degree of freedom system [31].For Example5, (a) comparing estimated and exact restoring forces, (b) displaying the estimation error of restoring force with m=51, and (c) displaying the estimation error of restoring force with m=201.
(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.



---

*Source: 841690-2013-11-05.xml*