# Research on Dynamic Load Identification Based on Explicit Wilson-θ and Improved Regularization Algorithm

Authors: Yuchuan Fan; Chunyu Zhao; Hongye Yu
Journal: Shock and Vibration
 (2019)

---

## Abstract

In the research of dynamic load identification, the method of obtaining kernel function matrix is usually rather cumbersome. To solve this problem, an explicit dynamic load identification algorithm based on the Wilson-θ (DLIAEW) method is proposed to easily obtain the kernel function matrix as long as the parameters of the system are known. To aim at the ill-posed problem in the inverse problem, this paper improves the Tikhonov regularization, proposes an improved regularization algorithm (IRA), and introduces the U-curve method to determine the regularization parameters. In the numeric simulation investigation of a four dofs vibrating system, effects of the sampling frequency and the noise level on the regularization parameters and the identification errors of impact and harmonic loads for the IRA are discussed in comparison with the Tikhonov regularization. Finally, the experiments of a cantilever beam excited by impact and harmonic loads are carried out to verify the advantages of the IRA.


---

## Body

## 1. Introduction

Load identification plays an important role in many engineering studies, such as reliability analysis, fault diagnosis, and health monitoring of mechanical power structures [1, 2]. Dynamic load identification methods include the frequency-domain method and the time-domain method [3]. Compared with the time-domain method, the study of dynamic load identification in the frequency domain starts earlier, and the theory is more mature. The frequency-domain method determines the dynamic force spectrum according to the relation between the transfer function matrix and the response spectrum of the system, or calculates the dynamic characteristics of the modal force in the frequency domain after the modal coordinate transformation [4]. The time-domain method is the inverse analysis based on the complex convolution relation between the load and the response, and the temporal history of dynamic loads is retrieved directly in the time domain. The time-domain method does not need Fourier transform, the result is intuitionistic, and the research of the time-domain method in recent years also has a great development.Dynamic load identification belongs to an ill-posedness of the inverse problem, which will lead to not useful solutions which cause large deviations from the exact solutions because of measured noise data and the randomness of structural parameters [5, 6]. In recent years, much effort on solving this ill-posed problem has been devoted to overcoming the effects of structural uncertainty and measurement noise and improving the accuracy of dynamic load identification. The effects of matrix ill-conditioning were overcome by using methods such as pseudoinversion of overdetermined matrixes, singular value rejection, singular value decomposition (SVD), and regularization techniques. Through both simulation and experiment on a flat rectangular plate, Thite and Thompson [7–9] proposed an assessment which was made of the success and failure of various strategies for dealing with the problems of ill-conditioning, in particular overdetermination and singular value rejection. Inoue et al. [10] utilized the SVD to locate the small singular values which were eliminated in computation of the frequency response function. Inoue et al. [11], Jacquelin et al. [12], and Adams and Doyle [13] described more systematic approaches of the SVD to solve the reconstruction problems of harmonic and nonharmonic forces.Using shape function to approximate dynamic load, kernel function response, and measured structure response, Liu et al. [14] established a time-domain dynamic Galerkin method (TDGM) to improve the accuracy of the identified dynamic load by taking shape function as the weighting function. Furthermore, they also proposed an efficient interpolation-based method to reduce ill-posedness by discretizing the load history into a series of time elements approximated through interpolation functions [15]. Qiao et al. [16–20] proposed the cubic B-spline collocation method and the sparse deconvolution method to identify impact loads. Through reconstructing the dynamic loads in the deterministic structure of the thin-walled cylindrical shell and airfoil structure, Wang et al. [21] developed a fast convergent iteration regularization method for identifying dynamic loads. The above works all focused on deterministic structures. However, practical structures have usually some random parameters. Hence, unknown load identified by virtue of a deterministic model describing a stochastic one would not be accurate. Expressing dynamic loads as the functions of time and random parameters in the time domain by Liu et al. [22, 23] transformed the dynamic models of uncertain structures into equivalent deterministic equations to identify unknown loads.In order to solve ill-posed problem for load identification, the choice of regularization parameters plays a key role in regularization, and there are many methods to select regularization parameters, such as the L-curve method, the generalized cross validation (GCV), the U-curve method, and the discrepancy principle [24–27]. To mitigate error propagation and ill-posed problem during the process of identification, Jia et al. [28, 29] proposed a weighted regularization approach based on the proper orthogonal decomposition (POD), in which the regularization parameter was selected by the GCV method.In this paper, an improved regulation algorithm is proposed to solve the ill-posed problem of dynamic load identification. In the next section, the explicit form of Wilson-θ is deduced, and a load identification algorithm based on the explicit Wilson-θ method is proposed. In Section 3, the IRA is proposed, and the U-curve method is introduced to select regularization parameters. The numeric simulations of a four dofs vibrating system and the experiments of a cantilever beam are conducted to verify the load identification effectiveness of the IRA in Sections 4 and 5, respectively. Finally, conclusions are given in Section 6.


## 2. Dynamic Load Identification Algorithm Based on Explicit Wilson-θ Method

The dynamic equation of a damped linear system with multiple dofs is(1)Mx¨+Cx˙+Kx=P,where M denotes the mass matrix, C denotes the damping matrix and K denotes the stiffness matrix. P is the load vector acting on the structure. x¨, x˙, and x denote the acceleration vectors, the velocity vectors, and the displacement vectors, respectively. Assuming that the structure is expressed as Rayleigh damping, it is expressed as(2)C=α1M+α2K,where α1 and α2 are the constants selected to achieve the desired damping ratio at two preselected periods/frequencies. The damping ratio for the nth mode of vibration is then given by [30](3)ξn=α12ωn+ωn2α2,where ωn corresponds to the circular frequency of the nth mode.In practical application, it is difficult to obtain the closed-form analytical solution of Equation (1) because of the complexity of structure and excitation. The dynamic equation is usually solved directly by a numerical integration method with the time-step method.The Wilson-θ method assumes that the acceleration varies linearly within the time interval t,t+θΔt θ≥1, as shown in Figure 1. In linear problems, the method is unconditionally stable for θ≥1.37, so θ=1.4 is usually employed [31].Figure 1
Model of the Wilson-θ method.Letδ be a time variable starting at time t, which is applicable to 0≤δ≤θΔt. According to the assumption of linear acceleration, the acceleration in this range can be expressed as follows:(4)x¨t+δ=x¨t+δθΔtx¨t+θΔt−x¨t.Integrating equation (4) once and twice, respectively, one can get(5)x˙t+δ=x˙t+x¨tδ+δ22θΔtx¨t+θΔt−x¨t,(6)xt+δ=xt+x˙tδ+12x¨tδ2+δ36θΔtx¨t+θΔt−x¨t.Insertingδ=θΔt into equations (5) and (6), we get(7)x˙t+θΔt=x˙t+θΔt2x¨t+θΔt+x¨t,(8)xt+θΔt=xt+θΔtx˙t+θ2Δt26x¨t+θΔt+2x¨t.According to equations (7) and (8), the acceleration and velocity of t+θΔt can be expressed by displacement:(9)x¨t+θΔt=6θ2Δt2xt+θΔt−xt−6θΔtx˙t−2x¨t,(10)x˙t+θΔt=3θΔtxt+θΔt−xt−2x˙t−θΔt2x¨t.Then, the dynamic equation att+θΔt can be expressed as follows:(11)Mx¨t+θΔt+Cx˙t+θΔt+Kxt+θΔt=Pt+θΔt,with(12)Pt+θΔt=Pt+θPt+Δt−Pt.Inserting equations (9), (10), and (12) into (11) yields(13)xt+θΔt=K^−1P^t+θΔt,with(14)K^=K+3θΔtC+6θ2Δt2M,(15)P^t+θΔt=Pt+θPt+Δt−Pt+M6θ2Δt2xt+6θΔtx˙t+2x¨t+C3θΔtxt+2x˙t+θΔt2x¨t.Solving equation (13) for xt+θΔt and then inserting xt+θΔt into equation (9) yields x¨t+θΔt. Substituting equation (9) into (4) and taking δ=Δt, we obtain(16)x¨t+Δt=6θ3Δt2xt+θΔt−xt−6θ2Δtx˙t+1−3θx¨t.Substituting equation (4) into (5) and (6) and taking δ=Δt, we get(17)x˙t+Δt=x˙t+Δt2x¨t+Δt+x¨t,(18)xt+Δt=xt+Δtx˙t+Δt26x¨t+Δt+2x¨t.Inserting equations (14) and (15) into (13) yields(19)xt+θΔt=K^−11−θPt+θK^−1Pt+Δt+6θ2Δt2MK^−1+3θΔtCK^−1xt+6θΔtMK^−1+2CK^−1x˙t+2MK^−1+θΔt2CK^−1x¨t.Substituting equation (19) into (16) and considering I=K^−1K^, we obtain(20)x¨t+Δt=C0Pt+C1Pt+Δt+Cdxt+Cvx˙t+Cax¨t,with(21)C0=6θ3Δt2K^−11−θ,C1=6θ2Δt2K^−1,Cd=−6θ3Δt2K^−1K,Cv=6θ2Δt6θ2Δt2MK^−1+2θΔtCK^−1−Ι,Ca=12θ3Δt2MK^−1+3θ2ΔtCK^−1+1−3θI.Inserting equation (20) into (17) yields(22)x˙t+Δt=B0Pt+B1Pt+Δt+Bdxt+Bvx˙t+Bax¨t,with(23)B0=Δt2C0,B1=Δt2C1,Bd=Δt2Cd,Bv=Δt2Cv+I,Ba=Δt2Ca+I.Inserting equation (20) into (18) yields(24)xt+Δt=A0Pt+A1Pt+Δt+Adxt+Avx˙t+Aax¨t,with(25)A0=Δt26C0,A1=Δt26C1,Ad=Δt26Cd+I,Av=Δt6Cv+IΔt,Aa=Δt26Ca+2I.Ifxi, x˙i, x¨i, and Pi represent the displacement, velocity, acceleration, and excitation in equations (20)∼(24) at time ti, the relationship of the system responses between the two adjoining time steps ti+1 and ti in the matrix form can be expressed as follows:(26)xi+1x˙i+1x¨i+1=A0A1B0B1C0C1PiPi+1+AdAvAaBdBvBaCdCvCaxix˙ix¨i,and the response at time ti can be also expressed as follows:(27)xix˙ix¨i=∑j=0i−1AdAvAaBdBvBaCdCvCajA0A1B0B1C0C1Pi−j−1Pi−j+AdAvAaBdBvBaCdCvCaix0x˙0x¨0,where i and j represent the power of the corresponding matrix, respectively. Equation (27) is a new explicit expression of the Wilson-θ method, and the responses of displacement, velocity, and acceleration in each time step can be solved simultaneously. Compared with the iterative algorithm, it is obvious that this explicit algorithm has more advantages. Letting Xi=xiT,xiT,xiTT and assuming that the system is zero initial responses, that is, X0=0 and P0=0, then equation (27) can be rewritten as follows:(28)Xi=∑j=0i−1DjEjPi−j−1Pi−j,i=1,2,⋯,n,with(29)Dj=AdAvAaBdBvBaCdCvCajA0B0C0,Ej=AdAvAaBdBvBaCdCvCajA1B1C1.Assuming thatyi denotes the measured responses of the system at time ti, we obtain(30)yi=LXi,where L is a mR×3n transformation matrix, in which mR is the number of the measured response, and n is the number of degrees of freedom. The element of the kth column for the lth row of the L matrix, which is corresponding to the lth measured response, is 1, and the other elements are 0. The value of k is(31)k=3×d−1+r,where d is corresponding to the dth degree of the measured response in the system and r is 1, 2, or 3, which represents that the measured response is displacement, velocity, and acceleration, respectively.If there aremF actual exciting forces in equation (1), that is, F=f1,f2,⋯,fmFT, which act on the l1th,l2th,⋯,lmFth masses, respectively, we obtain(32)Pi=RFi,where R is a n×mF transformation matrix, in which the element of the kth column for the lkth row of the R matrix is 1, and the other ones are 0, k=1,2,⋯,mF.Inserting equation (32) into (28) and left multiplying it by L, considering equation (30) and rearranging it, we obtain(33)yi=W0F1,i=1,∑j=1i−1Qi−j−1+Wi−jFj+W0Fi,i>1,with(34)Qj=LDjR,Wj=LEjR.Equation (33) can be written as a convolution form of matrices from time 1 to n:(35)Y=HF,with(36)Y=y1y2⋮yn,H=W000⋯0Q0+W1W00⋯0Q1+W2Q0+W1W0⋯0⋮⋮⋮⋱⋮Qn−2+Wn−1Qn−3+Wn−2Qn−4+Wn−3⋯W0,F=F1F2⋮Fn.For a deterministic system,H is a constant, and Y can be obtained by measuring the response of the system. In fact, the measurement data Y will be affected by the errors, such as measurement noise, truncation error, and principle error. Hence, the form of load identification with errors can be written as(37)Y+w=HF,where w denotes the error vector. The measurement data Y are disturbed by errors. At this time, if the condition number of the core matrix H is very large, the system will be seriously ill-posed. This means that the excitation F is very sensitive to even a small error disturbance of the measured data Y, so it is difficult to get an accurate excitation F. In order to overcome the influence of ill-posed problems, the regularization method can be used to determine F.


## 3. An Improved Regularization Algorithm

The regularization method is an effective means to solve ill-posed problems. The quality of the regularization method directly affects the accuracy of dynamic load identification results. Based on the Tikhonov regularization, an improved regularization method is proposed, and U-curve method is introduced to select the regularization parameter.


### 3.1. Tikhonov Regularization

In the load identification problem, the measurement data with noise,Yδ=Y+w, should satisfy(38)Yδ−Y≤δ.The purpose of Tikhonov regularization is to find the weighted minimum of the residual norm and solution norm in equation (39), so as to obtain the stable solution Fλ,δ of the ill-posed inverse problem [14]:(39)Fλ,δ=argminFHF−Yδ2+λF2,where HF−Yδ2 denotes the l2-norm of residuals, λ denotes the regularization parameter, and λF2 plays a role in regularizing ill-posed inverse. The norm term in equation (39) can be transformed into the matrix form:(40)JF,λ=HF−Yδ2+λF2=HF−YδTHF−Yδ+λFTF.The regular solution of load identificationFλ,δ can be obtained by finding the minimum of JF,λ:(41)∂JF,λ∂F=2HTHF−Yδ+2λF=0.The regular solutionFλ,δ should satisfy(42)HTH+λIFλ,δ=HTYδ.The explicit unique solution ofFλ,δ can be obtained from equation (42):(43)Fλ,δ=HTH+λI−1HTYδ.


### 3.2. Improved Regularization Algorithm

Inserting equation (32) into (26), left multiplying it by L, considering equation (30), we obtain(44)yi+1=Q0W0FiFi+1+yri,where yri denotes the responses of the initial condition for the system and is expressed as follows:(45)yri=LAdAvAaBdBvBaCdCvCaXi.InsertingFi+1=Fi+ΔFi into equation (44) and rearranging it, we obtain(46)W0ΔFi=yi+1−yri−yFi,with(47)yFi=W0+Q0Fi.Assuming that the measured response with noise isyi+1δ=yi+1+wi+1, and that ysi+1=yi+1−yFi−yri, the disturbance equation of (46) can be rewritten as follows:(48)W0ΔFi+ΔFni=ysi+1+wi+1,where ΔFni is the disturbance for the solution of equation (46).According to Reference [32], we obtain(49)ΔFniΔFi≤condW0wi+1ysi+1,where cond W0 is the condition number of the matrix W0 and expressed as follows:(50)condW0=W0⋅W0−1=maxW0zz⋅minW0zz−1.If condW0 is far greater than 1, effect of wi+1 on the calculated result is enhanced. On the other hand, the state variables Xi and force Fi need to be calculated from the measured responses y1δ,y2δ,⋯,yiδ during the process of load identification. Hence, ΔFi in equation (48) also consists of real increment, ΔFRi, of the exciting forces and pseudoincrement, ΔFPi, stemming from measurement noise and truncation error. The real increment of the exciting forces, ΔFRi, in interval of one sampling period is very small; that is, ΔFRi approaches to 0 with increase of the sampling frequency. According to equations (48) and (49), one can deduce that ΔFPi+ΔFni must be far greater than ΔFRi. Therefore, the least squares method of the ill-posed inverse problem Equation (37) can be expressed as follows:(51)minimizeHF−Yδ2,subject to1nΔF2≈0,where n denotes the sampling number. Classical penalty function methods [33] replace the constrained problem by an unconstrained problem of the form as follows:(52)minimizeHF−Yδ2+λ′1nΔF2.Assuming thatλ=λ′/n, we obtain(53)Fλ,δ=argminFHF−Yδ2+λΔF2.It is noticed that the solutionΔFi+ΔFni of equation (48) can also be obtained after solving equation (37) for F. Under the zero initial conditions, that is, X0=0 and F0=0, ΔFi+ΔFni can be expressed as follows:(54)ΔF=F1−F0F2−F1⋯Fn−Fn−1T=NF,with(55)N=1−11⋱⋱−11−11.Then, equation (53) can be transformed into(56)Fλ,δ=argminFHF−Yδ2+λNF2.Comparing equation (56) with equations (39) and (56) is considered to be an improved regularization algorithm (IRA). Similar to the Tikhonov regularization method, the norm term in equation (56) can be transformed into the matrix form:(57)JF,λ=HF−Yδ2+λF2=HF−YδTHF−Yδ+λNFTNF.The regular solution of load identificationFλ,δ can be obtained by finding the minimum of JF,λ:(58)∂JF,λ∂F=2HTHF−Yδ+2λNTNF=0.The regular solutionFλ,δ should satisfy(59)HTH+λNTNFλ,δ=HTYδ.The explicit unique solution ofFλ,δ can be obtained from the following equation:(60)Fλ,δ=HTH+λNTN−1HTYδ.


### 3.3. Determination of Regularization Parameters by U-Curve Method

The regularization method is actually a constrained least squares method. When the regularization parameterλ=0, equations (40) and (60) are transformed into the least squares method. The selection of regularization parameters is very important. In fact, ΔF∞ is always greater than 0. Hence, if the value of λ is too large, the recognition error will increase; if the value of λ is too small, it will not play the role of regularization.The principle of the U-curve method is similar to the L-curve method. The L-curve method [24] obtains the values of HF−Yδ and NF at different values of λ. Then, a curve with logHF−Yδ as the horizontal axis and logNF as the vertical axis is drawn. The shape of the curve is similar to the letter “L.” The optimum regularization parameter is the value of the position where the curvature of the curve is the largest (that is, the inflection point of the letter L). The U-curve method can be defined as follows [34]:(61)Uλ=1HF−Yδ2+1NF2,with(62)HF−Yδ2=∑i=1rλ4ψi2σi2+λ22,NF2=∑i=1rσi2ψi2σi2+λ22,where ψ=uTYδ, u is the left singular value matrix after the singular value decomposition of the kernel function matrix H. σ1≥σ2≥⋯≥σr>0 is the singular value of matrix H. The singular values of matrix N are all 1. Because the image of Uλ is similar to the letter “U,” it is known as the U-curve method. The purpose of the U-curve criterion is to select the parameters whose curvature reaches the local maximum of the left vertical part of the U-curve, where the value λ is the best regularization parameter. Because the U-curve method has less computation, it is superior to the L-curve method in calculating speed.


## 3.1. Tikhonov Regularization

In the load identification problem, the measurement data with noise,Yδ=Y+w, should satisfy(38)Yδ−Y≤δ.The purpose of Tikhonov regularization is to find the weighted minimum of the residual norm and solution norm in equation (39), so as to obtain the stable solution Fλ,δ of the ill-posed inverse problem [14]:(39)Fλ,δ=argminFHF−Yδ2+λF2,where HF−Yδ2 denotes the l2-norm of residuals, λ denotes the regularization parameter, and λF2 plays a role in regularizing ill-posed inverse. The norm term in equation (39) can be transformed into the matrix form:(40)JF,λ=HF−Yδ2+λF2=HF−YδTHF−Yδ+λFTF.The regular solution of load identificationFλ,δ can be obtained by finding the minimum of JF,λ:(41)∂JF,λ∂F=2HTHF−Yδ+2λF=0.The regular solutionFλ,δ should satisfy(42)HTH+λIFλ,δ=HTYδ.The explicit unique solution ofFλ,δ can be obtained from equation (42):(43)Fλ,δ=HTH+λI−1HTYδ.


## 3.2. Improved Regularization Algorithm

Inserting equation (32) into (26), left multiplying it by L, considering equation (30), we obtain(44)yi+1=Q0W0FiFi+1+yri,where yri denotes the responses of the initial condition for the system and is expressed as follows:(45)yri=LAdAvAaBdBvBaCdCvCaXi.InsertingFi+1=Fi+ΔFi into equation (44) and rearranging it, we obtain(46)W0ΔFi=yi+1−yri−yFi,with(47)yFi=W0+Q0Fi.Assuming that the measured response with noise isyi+1δ=yi+1+wi+1, and that ysi+1=yi+1−yFi−yri, the disturbance equation of (46) can be rewritten as follows:(48)W0ΔFi+ΔFni=ysi+1+wi+1,where ΔFni is the disturbance for the solution of equation (46).According to Reference [32], we obtain(49)ΔFniΔFi≤condW0wi+1ysi+1,where cond W0 is the condition number of the matrix W0 and expressed as follows:(50)condW0=W0⋅W0−1=maxW0zz⋅minW0zz−1.If condW0 is far greater than 1, effect of wi+1 on the calculated result is enhanced. On the other hand, the state variables Xi and force Fi need to be calculated from the measured responses y1δ,y2δ,⋯,yiδ during the process of load identification. Hence, ΔFi in equation (48) also consists of real increment, ΔFRi, of the exciting forces and pseudoincrement, ΔFPi, stemming from measurement noise and truncation error. The real increment of the exciting forces, ΔFRi, in interval of one sampling period is very small; that is, ΔFRi approaches to 0 with increase of the sampling frequency. According to equations (48) and (49), one can deduce that ΔFPi+ΔFni must be far greater than ΔFRi. Therefore, the least squares method of the ill-posed inverse problem Equation (37) can be expressed as follows:(51)minimizeHF−Yδ2,subject to1nΔF2≈0,where n denotes the sampling number. Classical penalty function methods [33] replace the constrained problem by an unconstrained problem of the form as follows:(52)minimizeHF−Yδ2+λ′1nΔF2.Assuming thatλ=λ′/n, we obtain(53)Fλ,δ=argminFHF−Yδ2+λΔF2.It is noticed that the solutionΔFi+ΔFni of equation (48) can also be obtained after solving equation (37) for F. Under the zero initial conditions, that is, X0=0 and F0=0, ΔFi+ΔFni can be expressed as follows:(54)ΔF=F1−F0F2−F1⋯Fn−Fn−1T=NF,with(55)N=1−11⋱⋱−11−11.Then, equation (53) can be transformed into(56)Fλ,δ=argminFHF−Yδ2+λNF2.Comparing equation (56) with equations (39) and (56) is considered to be an improved regularization algorithm (IRA). Similar to the Tikhonov regularization method, the norm term in equation (56) can be transformed into the matrix form:(57)JF,λ=HF−Yδ2+λF2=HF−YδTHF−Yδ+λNFTNF.The regular solution of load identificationFλ,δ can be obtained by finding the minimum of JF,λ:(58)∂JF,λ∂F=2HTHF−Yδ+2λNTNF=0.The regular solutionFλ,δ should satisfy(59)HTH+λNTNFλ,δ=HTYδ.The explicit unique solution ofFλ,δ can be obtained from the following equation:(60)Fλ,δ=HTH+λNTN−1HTYδ.


## 3.3. Determination of Regularization Parameters by U-Curve Method

The regularization method is actually a constrained least squares method. When the regularization parameterλ=0, equations (40) and (60) are transformed into the least squares method. The selection of regularization parameters is very important. In fact, ΔF∞ is always greater than 0. Hence, if the value of λ is too large, the recognition error will increase; if the value of λ is too small, it will not play the role of regularization.The principle of the U-curve method is similar to the L-curve method. The L-curve method [24] obtains the values of HF−Yδ and NF at different values of λ. Then, a curve with logHF−Yδ as the horizontal axis and logNF as the vertical axis is drawn. The shape of the curve is similar to the letter “L.” The optimum regularization parameter is the value of the position where the curvature of the curve is the largest (that is, the inflection point of the letter L). The U-curve method can be defined as follows [34]:(61)Uλ=1HF−Yδ2+1NF2,with(62)HF−Yδ2=∑i=1rλ4ψi2σi2+λ22,NF2=∑i=1rσi2ψi2σi2+λ22,where ψ=uTYδ, u is the left singular value matrix after the singular value decomposition of the kernel function matrix H. σ1≥σ2≥⋯≥σr>0 is the singular value of matrix H. The singular values of matrix N are all 1. Because the image of Uλ is similar to the letter “U,” it is known as the U-curve method. The purpose of the U-curve criterion is to select the parameters whose curvature reaches the local maximum of the left vertical part of the U-curve, where the value λ is the best regularization parameter. Because the U-curve method has less computation, it is superior to the L-curve method in calculating speed.


## 4. Numerical Example

In this section, a multi-degree-of-freedom vibrating system is used to verify the effectiveness of the proposed method in comparison with the Tikhonov regularization method. Through the numeric simulations, the impact and harmonic loads of the system are identified using the IRA and the Tikhonov regularization, respectively. The effects of the sampling frequency and noise level on the errors of load identification are discussed.The error between the recognition result and the real load is calculated through the following equation:(63)error1=Fid−FrealFreal×100%,error2=Fid−FrealFreal×100%,where error1 denotes the relative error, error2 denotes the relative accumulation error, Fid denotes the recognized load, Freal denotes the real force, ⋅ denotes the absolute value, and ⋅ denotes the l2-norm.Figure2 shows a four-degree-of-freedom system, in which masses are m1 = 4258 kg, m2 = 2258 kg, m3 = 2437 kg, and m4 = 1729 kg, respectively. The stiffness is k = 8.64 × 106 N/m. The damping is the Rayleigh damping, and the constants are α1 = 0.1 and α2 = 3.6157 × 10−4.Figure 2
Four-degree-of-freedom system.An impact load and the harmonic load are applied to the third degree of freedom of the system, respectively. The impact load is(64)ft=30e−2πp0t−0.82,where parameter p0 can adjust the frequency band of the impact load, and the peak value of the impact load appears at 0.8 s.Harmonic load is(65)ft=500sin10πt.In engineering, the measured responses of the system are inevitably disturbed by environmental noise. In order to simulate the measured data with noise, pseudorandom noise with uniform distribution is added to the displacement response. Hence, the displacement response is expressed as follows [35]:(66)Y˜=Y+e=Y+lnoise⋅stdy⋅2randn,1−1,where lnoise denotes the noise level, which is chosen as 10%, 20%, and 30%, respectively. std⋅ is a MATLAB script function that represents the standard deviation of vectors. The script function of MATLAB, rand (n, 1), can get a random vector of n×1, which represents the pseudorandom values in the standard uniform distribution on the interval of [0, 1].Using the Newmark-β algorithm, equation (1) can be numerically calculated to obtain displacement responses of the system. Time step is 0.001 s, and the calculating time is 3 s. Figures 3 and 4 show the displacement responses of m2 before and after adding 20% noise under the impact load and the harmonic load, respectively.Figure 3
Displacement response ofm2 for the system excited by the impact load before and after adding 20% noise.Figure 4
Displacement response ofm2 for the system excited by the harmonic load before and after adding 20% noise.Herein, the displacement responses with noise in Figures3 and 4 are used to identify the impact load and harmonic loads, respectively. The condition number of kernel function matrix H is 3.44e+23. In this case, matrix W0 is a scalar with a value of 1.095e − 12. According to equation (48), the error items, including the current measurement noise wi+1 and calculated errors in ysi+1, will be amplified by 9.13e11 times. These facts mean that the load identification problem of the system is seriously ill-posed. Using the U-curve method, the regularization parameter of the impact responses, λ, is determined to be 2e−12 and 3e−15 for the IRA and the Tikhonov regularization method, respectively; that of the harmonic response is 5e−13 and 6e−15 for the IRA and the Tikhonov regularization method, respectively. The condition numbers of the kernel function matrix regularized by the IRA and the Tikhonov regularization method are 1.88e+3 and 9.45e+3, respectively. Figure 5 shows comparisons of the indentified loads with the real ones, and the relative errors between the recognition result and the real load for the IRA and the Tikhonov regularization method are shown in Figure 6. As illustrated in Figures 5 and 6, the error of the identified load for the IRA is smaller than that for the Tikhonov regularization, except for the end segment of load identification.Figure 5
Loads identified by the responses with noise level of 20%: (a) impact load and (b) harmonic load.
(a)
(b)Figure 6
Relative errors identified by the responses with noise level of 20%: (a) impact load and (b) harmonic load.
(a)
(b)In order to draw a performance comparison between the IRA and the Tikhonov regularization, Tables1 and 2 list the regularization parameter and relative accumulation errors of identified loads for impact and harmonic loads under different sampling frequencies and different noise levels, respectively.Table 1
Regularization parameters and relative accumulation errors of the impact load.
Sampling frequency (Hz)
Noise level (%)
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
500
0
1e−15
1e−16
0.183671
0.230995
10
1e−13
1e−15
4.144848
8.688246
20
8e−13
1e−15
6.174314
12.484741
30
1e−12
2e−15
8.438775
19.757432
1000
0
1e−15
1e−16
0.064986
0.160739
10
1e−12
1e−15
2.557199
6.518704
20
2e−12
3e−15
5.923524
10.009010
30
2.7e−12
5e−15
9.163410
17.801819
2000
0
1e−15
1e−16
0.179296
2.679824
10
8e−13
1e−15
2.836595
3.899927
20
1e−12
5e−15
7.202332
8.112095
30
2e−12
7e−15
11.814816
12.270310Table 2
Regularization parameters and relative accumulation errors of harmonic load.
Sampling frequency (Hz)
Noise level (%)
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
500
0
1e−16
1e−17
0.951661
3.896035
10
1e−14
2e−15
3.087554
5.168342
20
5e−13
2e−15
5.093864
9.398055
30
1e−12
3e−15
7.489087
11.361723
1000
0
1e−16
1e−17
0.166175
0.20726
10
1e−14
2e−15
2.119125
4.458908
20
5e−13
6e−15
3.797509
7.188128
30
1e−12
3e−15
5.182875
10.441882
2000
0
1e−14
1e−16
0.099846
0.316657
10
1e−13
1e−15
2.02805
3.813177
20
5e−13
6e−15
2.582656
5.145480
30
4e−12
5e−15
3.340285
8.215345As shown in Tables1 and 2, the regularization parameter of the IRA is much bigger than that of the Tikhonov regularization method. Furthermore, the regularization parameters of the two algorithms all increase with increases of the sampling frequency and the noise level. But the difference of regularization parameters for the two algorithms also increases with increases of the sampling frequency and the noise level.The reason is that the single step increment of real exciting forces decreases with increase of the sampling frequency in comparison with the values of real exciting forces. Because eigenvalues of matrixN in equation (56) are the same as that of a unit matrix in equation (39), the condition number of kernel function matrix H regularized by the IRA must be much smaller than that regularized by the Tikhonov regularization method. Hence, relative accumulation errors of loads identified by the IRA are smaller than that by the Tikhonov regularization method. These facts demonstrate that noise resistance of the IRA is higher than that of the Tikhonov regularization and that the IRA is suitable to deal with complicated structures with huge sizes of mass and stiffness matrices as the same as the Tikhonov regularization method.


## 5. Experimental Verification

### 5.1. Experimental Setup

The experiments applied for load identification were conducted on an experimental system of vibration dynamics (model YE63251), as illustrated in Figure7. The system consists of an electric activator (model YE15400), an impact hammer (model LC-10A) including a force sensor (model CL-YD-303A) and an electric eddy current displacement sensor (model CWY-DO-502), and a cantilever beam. The parameters of the cantilever beam are listed in Table 3, and the first four natural frequencies of the cantilever beam are 15.9 Hz, 99.8 Hz, 279.5 Hz, and 547.7 Hz, respectively.Figure 7
Experimental system of cantilever beam: (a) impact load and (b) harmonic load.
(a)
(b)Table 3
Parameters of the cantilever beam.
Young’s modulus,E
200 GPa
Density,ρ
7840 kg/m3
Poisson’s ratio,μ
0.33
Length,L
0.64 m
Width,b
0.056 m
Height,h
0.008 m


### 5.2. Theoretical Model

In order to identify load of the beam using the corresponding response, the dynamic model is set up by virtue of FEM. The cantilever beam is assumed to be a plane one and divided into 18 segments, as shown in Figure8(b). As described by Hu and Wang [36], assuming that x-axis and y-axis represent the global coordinate system, the nodal displacements are grouped in the vector:(67)de=viθivjθjT,where v represents the translational dof and θ denotes the rotational dof.Figure 8
Structure diagram of cantilever beam: (a) planar beam element and (b) distribution of nodes for the cantilever beam.
(a)
(b)The element stiffness matrix (Ke) and the element mass matrix (Me) are expressed as follows:(68)Ke=EIzl3126l−126l6l4l2−6l2l2−12−6l12−6l6l2l2−6l4l2,Me=ρAJl42015622l54−13l22l4l13l−3l25413l156−22l−13l−3l2−22l4l2,where Iz denotes the moment of inertia, Iz=bh3/12; AJ denotes the cross-sectional area, AJ=bh; and l=L/18 denotes the length of the beam element.According to the finite element theory, the element stiffness matrixKe and the element mass matrix Me are integrated into the total stiffness matrix K and the total mass matrix M, respectively [36]. Then, the dynamic equations of the system can be expressed as follows:(69)Md¨+α1M+α2Kd˙+Kd=0.Using matricesM and K, the natural frequencies of the beam were determined. The first five natural frequencies are listed in the second row of Table 4.Table 4
Comparison of natural frequencies for technical manual, theoretical, and measured values (Hz).
Mode
1st
2nd
3rd
4th
Technical manual
15.9
99.8
297.5
547.7
Theoretical values
15.7
99.4
296.9
547.1
Measured values
15.6
98.9
—
—In order to determine the parameters of the finite element model, a hammer impact test was conducted to record impact load and displacement response, as shown in Figure7(a). The hammer excited the beam at node 10, and the displacement sensor was located at node 12. Figure 9 shows the response of the beam for the hammer impact and its FFT results. The hammer impact excited only the first two natural frequencies, 15.6 Hz and 98.9 Hz, as shown in Figure 9(b). These values were slightly smaller that their theoretical ones, as shown in Table 4. The fact stems from the contact stiffness of clamped end of the cantilever beam. As shown in Figure 9(a), the amplitude enclosure line of the response was approximately an exponential decay. This fact demonstrates that damping of the system is Rayleigh damping. According to the experimental results in [37], the hammer impact excites the free vibration frequencies of the system, which can be expressed as follows:(70)yt=A0e−ξωntsinωdt+γ,with(71)ωd=ωn1−ξ2.Figure 9
Impact response of the cantilever beam and its FFT results: (a) impact response and (b) FFT results.
(a)
(b)In comparison with the amplitude of component 15.6 Hz, that of component 98.9 Hz is very small and can be neglected. So, the relationship between the amplitude ofyt and time t can be expressed as follows:(72)Ampyt=A0e−ξ1ωn1t,orlnAmpyt=a0−b0t,where a0=lnA0 and b0=ξ1ωn1.The amplitude data were extracted from the displacement response in Figure10(a) in the range of 6 s to 38 s. Then, the data of lnAmpyt were linearly fitted using the least square method, and the results were shown in Figure 10. The coefficients of the linear fitting were a0=6.04 and b0=0.1117, that is, ξ1ωn1=0.1117. Using equation (71), the natural frequency and the critical damping ratio of the first mode for the cantilever beam were determined to be 15.6 Hz and 0.00114, respectively. Taking the natural frequency of the second mode as 99.4 Hz, the critical damping ratio of the second mode is determined to be ξ2=0.1002 using equation (71).Figure 10
Curve fittings of amplitude for the impact response: (a) power function curve fitting and (b) linear fitting.
(a)
(b)Using equation (3), we get linear equations as follows:(73)α162.4π+15.6πα2=0.00114,α1397.6π+99.4πα2=0.1002.Solving equation (73) for α1 and α2, we obtain(74)α1=−2.9315,α2=3.2839×10−4.In order to ensure the fidelity between the finite element model and the actual specimen, a harmonic exciting test was also conducted to record exciting force and the displacement response, as shown in Figure7(b). The electric activator excited the beam at node 10, and the displacement sensor was located at node 12. Using the Newmark-β algorithm, Equation (69) was simulated to determine the displacement responses excited by the impact load and the harmonic load, respectively. Figures 11(a) and 11(b) show comparisons of the simulation result with the experimental result for the impact and harmonic loads, which reveal quite well agreement, respectively.Figure 11
Comparisons of the simulation result with the experimental result: (a) impact load and (b) harmonic load.
(a)
(b)


### 5.3. Load Identifications

Through the data of displacement responses in Figure11, the impact and harmonic loads of the system were identified using the IRA and the Tikhonov regularization method, respectively. The identified results of impact load are shown in Figure 12, and that of harmonic load and its FFT result are shown in Figure 13. Furthermore, the measured point of displacement response stayed at the point of node 12, tests of excited point at the points of nodes 4, 6, and 15 for impact and harmonic loads were conducted to record the exciting force and the corresponding displacement response, respectively. Then, the excited loads were identified using the IRA and the Tikhonov regularization method, respectively.Figure 12
Identified results of the impact load.Figure 13
Identified result of the harmonic loads and its frequency spectrum.
(a)
(b)
(c)Tables5 and 6 show comparisons of the IRA with the Tikhonov regularization method for the regularization parameters and the relative accumulation errors of impact and harmonic loads, respectively. As the same as the simulation results in Tables 1 and 2, the regularization parameter of the IRA is much bigger than that of the Tikhonov regularization method and the relative accumulation errors of loads identified by the virtue of the IRA are smaller than that by the virtue of the Tikhonov regularization. All the above facts demonstrate that the performance of the IRA proposed in this paper is better than that of the Tikhonov regularization method for load identification.Table 5
Comparisons of the IRA with the Tikhonov regularization for the identified impact load errors under differential exciting points.
Excitation points
Condition number
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
6
4.7258e6
6.2546e−9
3.6244e−11
4.651456
8.558625
8
4.9454e6
3.2258e−9
2.6716e−11
3.869557
7.776456
10
5.5213e6
4.5743e− 9
1.5416e−11
3.568851
7.216584
15
5.8876e6
6.2538e− 9
5.2264e−11
4.655242
8.588946Table 6
Comparisons of the IRA with the Tikhonov regularization for the identified harmonic load errors under differential excited points.
Excitation points
Condition number
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
4
4.5079e6
7.4715e−11
2.8212e−13
3.298616
4.298616
10
5.1510e6
2.9021e−9
1.2324e−11
3.414175
4.625209
12
5.5405e6
3.5235e−9
1.4592e−11
3.239187
4.252354
15
6.1644e6
8.8048e−9
3.7205e−11
2.348616
5.474881


## 5.1. Experimental Setup

The experiments applied for load identification were conducted on an experimental system of vibration dynamics (model YE63251), as illustrated in Figure7. The system consists of an electric activator (model YE15400), an impact hammer (model LC-10A) including a force sensor (model CL-YD-303A) and an electric eddy current displacement sensor (model CWY-DO-502), and a cantilever beam. The parameters of the cantilever beam are listed in Table 3, and the first four natural frequencies of the cantilever beam are 15.9 Hz, 99.8 Hz, 279.5 Hz, and 547.7 Hz, respectively.Figure 7
Experimental system of cantilever beam: (a) impact load and (b) harmonic load.
(a)
(b)Table 3
Parameters of the cantilever beam.
Young’s modulus,E
200 GPa
Density,ρ
7840 kg/m3
Poisson’s ratio,μ
0.33
Length,L
0.64 m
Width,b
0.056 m
Height,h
0.008 m


## 5.2. Theoretical Model

In order to identify load of the beam using the corresponding response, the dynamic model is set up by virtue of FEM. The cantilever beam is assumed to be a plane one and divided into 18 segments, as shown in Figure8(b). As described by Hu and Wang [36], assuming that x-axis and y-axis represent the global coordinate system, the nodal displacements are grouped in the vector:(67)de=viθivjθjT,where v represents the translational dof and θ denotes the rotational dof.Figure 8
Structure diagram of cantilever beam: (a) planar beam element and (b) distribution of nodes for the cantilever beam.
(a)
(b)The element stiffness matrix (Ke) and the element mass matrix (Me) are expressed as follows:(68)Ke=EIzl3126l−126l6l4l2−6l2l2−12−6l12−6l6l2l2−6l4l2,Me=ρAJl42015622l54−13l22l4l13l−3l25413l156−22l−13l−3l2−22l4l2,where Iz denotes the moment of inertia, Iz=bh3/12; AJ denotes the cross-sectional area, AJ=bh; and l=L/18 denotes the length of the beam element.According to the finite element theory, the element stiffness matrixKe and the element mass matrix Me are integrated into the total stiffness matrix K and the total mass matrix M, respectively [36]. Then, the dynamic equations of the system can be expressed as follows:(69)Md¨+α1M+α2Kd˙+Kd=0.Using matricesM and K, the natural frequencies of the beam were determined. The first five natural frequencies are listed in the second row of Table 4.Table 4
Comparison of natural frequencies for technical manual, theoretical, and measured values (Hz).
Mode
1st
2nd
3rd
4th
Technical manual
15.9
99.8
297.5
547.7
Theoretical values
15.7
99.4
296.9
547.1
Measured values
15.6
98.9
—
—In order to determine the parameters of the finite element model, a hammer impact test was conducted to record impact load and displacement response, as shown in Figure7(a). The hammer excited the beam at node 10, and the displacement sensor was located at node 12. Figure 9 shows the response of the beam for the hammer impact and its FFT results. The hammer impact excited only the first two natural frequencies, 15.6 Hz and 98.9 Hz, as shown in Figure 9(b). These values were slightly smaller that their theoretical ones, as shown in Table 4. The fact stems from the contact stiffness of clamped end of the cantilever beam. As shown in Figure 9(a), the amplitude enclosure line of the response was approximately an exponential decay. This fact demonstrates that damping of the system is Rayleigh damping. According to the experimental results in [37], the hammer impact excites the free vibration frequencies of the system, which can be expressed as follows:(70)yt=A0e−ξωntsinωdt+γ,with(71)ωd=ωn1−ξ2.Figure 9
Impact response of the cantilever beam and its FFT results: (a) impact response and (b) FFT results.
(a)
(b)In comparison with the amplitude of component 15.6 Hz, that of component 98.9 Hz is very small and can be neglected. So, the relationship between the amplitude ofyt and time t can be expressed as follows:(72)Ampyt=A0e−ξ1ωn1t,orlnAmpyt=a0−b0t,where a0=lnA0 and b0=ξ1ωn1.The amplitude data were extracted from the displacement response in Figure10(a) in the range of 6 s to 38 s. Then, the data of lnAmpyt were linearly fitted using the least square method, and the results were shown in Figure 10. The coefficients of the linear fitting were a0=6.04 and b0=0.1117, that is, ξ1ωn1=0.1117. Using equation (71), the natural frequency and the critical damping ratio of the first mode for the cantilever beam were determined to be 15.6 Hz and 0.00114, respectively. Taking the natural frequency of the second mode as 99.4 Hz, the critical damping ratio of the second mode is determined to be ξ2=0.1002 using equation (71).Figure 10
Curve fittings of amplitude for the impact response: (a) power function curve fitting and (b) linear fitting.
(a)
(b)Using equation (3), we get linear equations as follows:(73)α162.4π+15.6πα2=0.00114,α1397.6π+99.4πα2=0.1002.Solving equation (73) for α1 and α2, we obtain(74)α1=−2.9315,α2=3.2839×10−4.In order to ensure the fidelity between the finite element model and the actual specimen, a harmonic exciting test was also conducted to record exciting force and the displacement response, as shown in Figure7(b). The electric activator excited the beam at node 10, and the displacement sensor was located at node 12. Using the Newmark-β algorithm, Equation (69) was simulated to determine the displacement responses excited by the impact load and the harmonic load, respectively. Figures 11(a) and 11(b) show comparisons of the simulation result with the experimental result for the impact and harmonic loads, which reveal quite well agreement, respectively.Figure 11
Comparisons of the simulation result with the experimental result: (a) impact load and (b) harmonic load.
(a)
(b)


## 5.3. Load Identifications

Through the data of displacement responses in Figure11, the impact and harmonic loads of the system were identified using the IRA and the Tikhonov regularization method, respectively. The identified results of impact load are shown in Figure 12, and that of harmonic load and its FFT result are shown in Figure 13. Furthermore, the measured point of displacement response stayed at the point of node 12, tests of excited point at the points of nodes 4, 6, and 15 for impact and harmonic loads were conducted to record the exciting force and the corresponding displacement response, respectively. Then, the excited loads were identified using the IRA and the Tikhonov regularization method, respectively.Figure 12
Identified results of the impact load.Figure 13
Identified result of the harmonic loads and its frequency spectrum.
(a)
(b)
(c)Tables5 and 6 show comparisons of the IRA with the Tikhonov regularization method for the regularization parameters and the relative accumulation errors of impact and harmonic loads, respectively. As the same as the simulation results in Tables 1 and 2, the regularization parameter of the IRA is much bigger than that of the Tikhonov regularization method and the relative accumulation errors of loads identified by the virtue of the IRA are smaller than that by the virtue of the Tikhonov regularization. All the above facts demonstrate that the performance of the IRA proposed in this paper is better than that of the Tikhonov regularization method for load identification.Table 5
Comparisons of the IRA with the Tikhonov regularization for the identified impact load errors under differential exciting points.
Excitation points
Condition number
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
6
4.7258e6
6.2546e−9
3.6244e−11
4.651456
8.558625
8
4.9454e6
3.2258e−9
2.6716e−11
3.869557
7.776456
10
5.5213e6
4.5743e− 9
1.5416e−11
3.568851
7.216584
15
5.8876e6
6.2538e− 9
5.2264e−11
4.655242
8.588946Table 6
Comparisons of the IRA with the Tikhonov regularization for the identified harmonic load errors under differential excited points.
Excitation points
Condition number
Regularization parameter
Error (%)
IRA
Tikhonov
IRA
Tikhonov
4
4.5079e6
7.4715e−11
2.8212e−13
3.298616
4.298616
10
5.1510e6
2.9021e−9
1.2324e−11
3.414175
4.625209
12
5.5405e6
3.5235e−9
1.4592e−11
3.239187
4.252354
15
6.1644e6
8.8048e−9
3.7205e−11
2.348616
5.474881


## 6. Conclusions

This paper gives a general explicit form of the Wilson-θ algorithm to identify dynamic loads. Then, an improved regularization algorithm is proposed to solve the ill-posed problem of dynamic load identification instead of the Tikhonov regularization method. The main novelty of the IRA is that it replaces l2-norm of the exciting force vector with that of the increment vector for the exciting force. Compared with the Tikhonov regularization method through the computer simulations and the experiments for identifications of impact and harmonic loads, the regularization parameter of the IRA is much bigger than that of the Tikhonov regularization method. Furthermore, the regularization parameters of the two algorithms all increase with increases of the sampling frequency and the noise level. But the difference of regularization parameters for the two algorithms also increases with increases of the sampling frequency and the noise level. Therefore, the performance of the IRA is better than that of the Tikhonov regularization method for load identification. Because all the eigenvalues of the transformation matrix between the exciting force vector and its increment vector are the same as that of an unit matrix, the IRA is suitable to deal with complicated structures with huge sizes of mass and stiffness matrices as the same as the Tikhonov regularization method.



---

*Source: 8756546-2019-03-25.xml*