Shock and Vibration

Volume 2016 (2016), Article ID 7328969, 16 pages

http://dx.doi.org/10.1155/2016/7328969

## A New Method for Optimal Regularization Parameter Determination in the Inverse Problem of Load Identification

^{1}Department of Astronautic Science and Mechanics, Harbin Institute of Technology (HIT), P.O. Box 304, No. 92 West Dazhi Street, Harbin 150001, China^{2}School of Mathematics and Statistics, Northeast Petroleum University, Daqing 163318, China

Received 9 December 2015; Accepted 1 February 2016

Academic Editor: Nerio Tullini

Copyright © 2016 Wei Gao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

According to the regularization method in the inverse problem of load identification, a new method for determining the optimal regularization parameter is proposed. Firstly, quotient function (QF) is defined by utilizing the regularization parameter as a variable based on the least squares solution of the minimization problem. Secondly, the quotient function method (QFM) is proposed to select the optimal regularization parameter based on the quadratic programming theory. For employing the QFM, the characteristics of the values of QF with respect to the different regularization parameters are taken into consideration. Finally, numerical and experimental examples are utilized to validate the performance of the QFM. Furthermore, the Generalized Cross-Validation (GCV) method and the* L*-curve method are taken as the comparison methods. The results indicate that the proposed QFM is adaptive to different measuring points, noise levels, and types of dynamic load.

#### 1. Introduction

Load and system identification are two kinds of problems of the inverse problem theory [1–3]. System and load identification play an important role in the fields of the stability analysis of the structure, structural health monitoring, and so forth [4, 5]. It is hard or impossible to measure dynamic load directly under operating condition for technical or economic reasons such as dynamic loads acting on high buildings and impact loads acting on aircraft. Load identification technique can be mainly divided into frequency-domain and time-domain methods [6, 7]. The frequency-domain method has limitations in the nonstationary and transient cases. Therefore, time-domain method which is based on the kinetic equation and the Duhamel integral equation has become the main focus [8–11]. A discrete linear system equation can be obtained as with the discretization of Duhamel integral equation [12], where denotes the system matrix. The Singular Value Decomposition (SVD) of the system matrix is necessary to solve the discrete linear system equations. However, SVD may cause ill-conditioned problems as the minimum singular value of the system matrix approaches zero. Obviously, it is necessary to solve the ill-conditioned problems resulting from SVD. The Tikhonov regularization method was proposed to overcome the ill-conditioned problems of load identification [13–15]. Subsequently, a lot of research literatures have focused on load identification problems based on the regularization method. Consequently, a series of load identification theories were introduced to satisfy various engineering problems in practice. Mao et al. proposed a method for load identification based on Markov parameters precise computation [16]. Liu et al. proposed a time-domain discretization load identification method which was based on regularization technology [17]. Sanchez and Benaroya reviewed some in-depth researches of regularization methods which had been applied on load identification successfully [18]. However, they did not discuss the determination of the optimal regularization parameter. Tikhonov regularization method of load identification can be divided into three parts: the establishment of the discrete linear system equation, the solution of Tikhonov regularization method, and the determination of the optimal regularization parameter. One principal direction of the Tikhonov regularization method for load identification is the determination of the optimal regularization parameter. Different methods, such as ridge trace method [19], the quasi optimal criterion [20], Generalized Cross-Validation (GCV) [21], Morozov discrepancy principle [22], heuristic rules [23], and* L*-curve criterion [24], can be utilized to determine the optimal regularization parameter.

In general, methods for determining the optimal regularization parameter can be divided into prior and posterior strategies. Compared with prior strategy, the posterior strategy does not need to find out the noise level of structural response in advance. The GCV method and the* L*-curve method are the two commonly used posterior strategy methods to select the optimal regularization parameter in engineering [16, 17, 25–28]. Although the convergence of the GCV criterion has been proved theoretically, the curve of the GCV function used to choose the optimal regularization parameter is sometimes too flat to find the minimum of the GCV function [29]. The* L*-curve method needs to ensure the approximate range of the optimal regularization parameter in advance [30]. Fortunately, the proposed QFM does not need to ensure the approximate range of the optimal regularization parameter in advance. The optimal regularization parameter can be determined easily and effectively.

This paper focuses on the regularization method in the inverse problem of load identification. The rest of the paper is organized as follows. In Section 2, based on appropriate improvement of the shape function method of moving least square fitting [17], the discrete linear system equation model for load identification is constructed. In Section 3, the Tikhonov regularization method for solving ill-conditioned problems is introduced. In Section 4, the GCV method and the* L*-curve method for determining the optimal regularization parameter are introduced. Meanwhile, based on quadratic programming theory, the QFM for determining the optimal regularization parameter is proposed. In Section 5, numerical simulation examples are used to illustrate the excellent adaptability of the QFM to different measuring points, noise levels, and types of dynamic load. Then the comparative results with the GCV method and the* L*-curve method are given. The proposed method is implemented to obtain the identification results of dynamic load acting on a cantilever beam structure in Section 6. Finally, some useful conclusions are drawn in Section 7.

#### 2. Discrete Linear System Equation Model

Based on the discretization of Duhamel integralthe corresponding discrete linear system equation can be obtained, where is a structural response which can be displacement, velocity, acceleration, and so forth. is the corresponding Green kernel function which is the impulse response of the structure. is the dynamic load acting on the structure which needs to be identified. Currently, many methods treat the load as a constant in each time unit. The real load can be expressed by the step function approximately. However, it is unable to represent the tendency of the real load at each time unit effectively. Moreover, there is almost no relativity between the approximate loads at the adjacent time units. According to the shape function method of moving least square fitting [17], the way for constructing the optimal approximate load was improved. The whole time domain is separated into many time units with length of . We suppose there are time units in all and are the loading points to be identified. Suppose is a local support domain in the whole time domain. The local support domain contains time units and loading points. Without losing generality, suppose these loading points in the local support domain are and are the basis functions. The approximate load expression can be written aswhere is the basis function vector, and is the combination coefficient vector. The residual between approximate load and real load in isIf we want to obtain the expression of the optimal approximate load, the number of loading points should be equal to or larger than that of basis functions . The necessary condition of the least residual is that the partial derivative of the combination coefficients in (3) is equal to zero. ThenthenLetthenwhere . Substituting (7) into (2) giveswhere is the shape function vector. Equation (8) is the optimal approximation of the real load in the local domain under some meaning of weighted least square. Then the approximate load in the first time unit of (8) is the optimal approximation of the real load between the loading points and . The corresponding shape function vector iswhere denotes zero vector with the same length as the shape function vector. When the length of is fixed, the shape function vector in (9) only depends on the basis function , and the shape function can be used in all time units with the same form. Thus, the optimal approximation of the real load in whole time domain can be obtained gradually via moving this supported domain .

If we use the approximate load in the first time unit of (8) to construct the optimal approximation of the real load in whole time domain, the shape function can be defined as follows:The expression of optimal approximate load in the whole time domain iswhere denotes the approximate load in time unit and . However, the loading points on the right side of do not exist in practice. Without losing generality, we suppose . ThenSubstituting (11) into (1) givesEquation (13) can be rearranged aswhere denotes the component of . From (14), the finite element computational amount to solve the response resulting from shape function is only times. Suppose that, through finite element computing, the structural responses caused by each shape function load areSimilar to the method of the assembly of stiffness matrixes in finite element method, the whole response matrix of shape functions can be formed by the means of assembling the responses of the shape function loads in all time units. The discretization of (14) givesWith the consideration of (12) and (14), the vector of all should beSimilarly, the discretization of (14) gives the shape function response matrixes asLetwhere denotes the column of . Then we have . From (17) and (19), letthen is the discrete linear system equation when we use the approximate load in the first time unit of (8) to construct the optimal approximation of the real load in whole time domain. The improved method to construct the discrete linear system equation can be established as follows.

Without losing generality, we suppose . Then we can obtain the shape function vector corresponding to the middle and the last time unit of in the local domain asFrom (9) to (11), the improved expression of optimal approximate load can be proposed as where denotes the optimal approximate load of the time unit. Meanwhile, , , , . Without losing generality, we suppose . Substituting (23) into (1) givesSimilar to the construction of , the shape function response matrix of and can be obtained with the discretization of (24) aswhere denote the column of . Letwhere denotes the component of and denotes the component of . Letwhere denotes zero matrix. Let . According to (16) and (21), the discrete linear system equation model can be obtained asEquation (28) is the model of discrete linear system equation based on single-input and single-output (SISO) system. For the multi-input and multioutput (MIMO) system, we suppose and are the two loads acting on different locations. The corresponding measured structural responses in two different locations are and . Therefore, according to (1), there are and . is the corresponding Green kernel function which is the impulse response of the structure. It is similar to the deduction process of (28); there are and , where can be obtained by the discretization of . Matrix denotes the system matrix. Then the discrete linear system equation model of MIMO system can be constructed as

#### 3. Tikhonov Regularization Method for Ill-Conditioned Problem

The dynamic load identification problem is ill-conditioned, because the system matrix in (28) is always ill-conditioned. Furthermore, there exists noise in the measured structure response. Therefore, it is almost impossible to get the accurate solution. The approximate numerical solution that satisfies the needs in practice can be obtained only. Now the Tikhonov regularization method is one of the most commonly used methods to solve the ill-conditioned problem. The main idea is to study the minimization problem aswhere is a constant and denotes Euler norm. The equivalent form of objective function in (30) iswhere is identity matrix. Using the idea that the gradient of objective function equals zero, the expression of the least square solution of (30) can be obtained asFrom SVD and [30], (32) can be rearranged aswhere is the real response data, is the corresponding noise data, and are all the singular values of matrix . Sometimes the system matrix is seriously ill-conditioned. Therefore, the value corresponding to the real response of the right side of (33) fluctuates heavily with the decrease of the singular value. In addition, its stability is also poor. Serious effects on the identified load will be caused by the measured noise , which will make the error uncontrollable between the identified load and the real load. However, the regularization operator plays an important role in (33). The identified load that posed both the stability of numerical inversion and the controllability of error with real load can be obtained by choosing a reasonable regularization parameter . Therefore, the determination of the optimal regularization parameter is the key problem of regularization method of load identification.

#### 4. Methods for Determining Optimal Regularization Parameter

##### 4.1. Generalized Cross-Validation (GCV) [21]

According to the discrete linear system equation in (28), the GCV function is defined aswhere denotes the trace of matrix. The determination of optimal regularization parameter will satisfy the condition

The image of GCV function is sometimes too flat to find the minimum of the GCV function [29]. Because the values of GCV function almost do not change near the minimum value, this brings us more trouble to determine its minimum and the optimal regularization parameter. In addition, the GCV method occasionally fails to determine the optimal regularization parameter [31]. The value of optimal regularization parameter obtained by GCV method is sometimes too large, which will result in losing too much information of the real response.

##### 4.2.
*L*-Curve Criterion [24]

The* L*-curve criterion utilizes log-log as the scale to describe the comparison between the norm of regularization solution and residual and adopts the biggest curvature of* L*-curve under the log-log scale as the criterion for determining optimal regularization parameter. Let and . The* L*-curve curvature isThe condition for determining optimal regularization parameter satisfies

The* L*-curve method needs to ensure the approximate range of optimal regularization parameter in advance [30]. Furthermore, there is no effective method for* L*-curve method to confirm optimal regularization parameter in its searching range because it is relatively small. We need several attempts to find ideal* L*-curve [31]. However, the larger searching range will result in the overlarge calculation burden or the unsatisfactory image of* L*-curve.

##### 4.3. Quotient Function Method (QFM)

The minimization problem in (30) can be solved by the quadratic programming algorithm [32–34]. When the value of regularization parameter is large, the matrix in (31) will be well-posed. Although matrixes and are both ill-conditioned, the matrix is well-posed when the value of regularization parameter is large enough. As a result, a good least square solution that satisfies (32) can be obtained asthen (38) can be rearranged asFrom (33), can be regarded as function of regularization parameter . Then the two sides of (39) can be regarded as the function of parameter . Then (39) can be rearranged asWith the value of regularization parameter decreasing, matrix will transform from being well-posed to being weakly ill-conditioned. Thus the least square solution that is solved through the quadratic programming algorithm is not fairly accurate. Then (40) becomesIf the value of regularization parameter decreases more, matrix will transform to being strongly ill-conditioned. The accuracy of the least square solution will become even worse. Then (41) becomesTherefore, quotient function is defined asWe hope to find out a parameter satisfying the following conditions:(i)if , then ;(ii)if , then .Then the parameter will be selected as the optimal regularization parameter. According to the theoretical analysis of QFM, the optimal regularization parameter can be determined effectively through the values of the QF corresponding to different regularization parameters.

##### 4.4. The Algorithm Realization of the QFM

(1)Get the largest singular value by performing SVD to the system matrix . The initial reference parameter vector can be defined as .(2)Calculate values of the QF corresponding to every component in the vector . The vector that is made of the values of QF can be obtained as .(3)Set a threshold that equals 1 approximately. All values of QF in vector are checked. If the distance between any component and 1 is less than or equal to , then the reference parameter vector can be renewed as .(4)Repeat steps () and () until the distance between a certain component of vector and 1 is more than (suppose the component is the component of ). Now we suppose the reference parameter vector is .(5)Stop the searching process with the optimal regularization parameter being the component of ; that is, .

The algorithm realization of the QFM can be explained as follows: the optimal regularization parameter is smaller than the largest singular value . Thus, the largest component of the initial reference vector is chosen as the largest singular value , which can guarantee the value of the optimal regularization parameter within the searching range. According to the different demands for accuracy in practice application, the threshold can be chosen as , , or . If there is no demand, we choose threshold as . According to the threshold , the searching process will be stopped if the distance between any component of vector and 1 is more than . Then we can determine the optimal regularization parameter from the conditions (i) and (ii) in Section 4.3.

#### 5. Numerical Simulation

A finite element model of cantilever beam is built considering the real beam structure to be introduced in Section 6, as shown in Figure 20. From left to right are nodes 1 to 10, as shown in Figure 1. Take , , and as the external load applied on node 8. Amplitude modulation (AM) load and sine load are acting on node 8 for 0.3 s. Sine load is acting on node 8 for 0.005 s to simulate impact load.