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.
In the numerical examples, is selected as the basis function vector to construct the system matrix in (28). The number of loading points . We can solve (28) by using the Tikhonov regularization method and determine the optimal regularization parameter through the QFM and the GCV and the L-curve methods. Three different types of dynamic load are identified at different noise levels. The acceleration response for load identification is assumed to suffer from the measured noise due to environmental effects. It is modeled with white noise that is added to the numerical calculated responses to simulate the polluted measurements aswhere is the acceleration response for load identification; is the numerical calculated acceleration response; is the noise level; is a vector of white noise with zero mean and unit standard deviation; denotes the standard deviation. The relative error (RE) and correlation coefficient (CC) between the identified load and the real load can be defined aswhere denotes the component of , denotes the real load, and denotes the average.
5.1. SISO System Load Identification in Numerical Simulation
5.1.1. Identification of the AM Load
The sampling frequency is 1000 Hz. Acceleration response of node 6 is selected for AM load identification at noise levels 1%, 3%, and 5%. The threshold is chosen as . At noise level 1%, the objective function of the minimization problem in (30) is constructed according to (28). The optimal regularization parameter is determined by QFM. The curve of QF is shown in Figure 2. Figure 2 indicates that when , the values of QF are far away from 1. When , almost all the values of QF are equal to 1, which satisfies the needs of QFM algorithm. Therefore, we choose the range of that satisfies the condition when QF are almost equal to 1. According to conditions (i) and (ii), the optimal regularization parameter is selected as the smallest value in the range of . As a result, the value of optimal regularization parameter determined by QFM is .
Similarly, the regularization parameter is based on the GCV method as shown in Figure 3. The result of L-curve method is as shown in Figure 4. The red stars in Figures 2, 3, and 4 denote the corresponding location of the optimal regularization parameter. The corresponding identified loads of three methods are shown in Figure 5. The relative results of identified load at different noise levels are shown in Table 1. Figure 5 shows that the identified loads of the three methods fit the value and reflect the variation trend of the real load quite well. The results of Table 1 illustrate that all these three methods can identify AM load effectively and obtain the stable solution with high accuracy.
5.1.2. Identification of the Impact Load
The sampling frequency is 5000 Hz. Acceleration response during 0–0.03 s of node 6 is selected for impact load identification at noise levels 1%, 3%, and 5%. The threshold is chosen as . The main characteristic parameter of impact load is the maximum amplitude. Thus, the RE of identified impact load is the error of the maximum amplitude between identified load and real load. At noise level 1%, the curves of QF and GCV and L-curve are shown in Figures 6, 7, and 8; the corresponding loads identified by these three methods are shown as Figure 9. At different noise levels, the relative results of identified impact load (IIL) are illustrated in Table 2.
Figures 6 and 9 indicate that the QFM can determine the optimal regularization parameter effectively. The corresponding identified load can well reflect the time history of impact load. The maximum amplitudes of identified load are nearly equal to real impact load. Table 2 shows that the three methods can obtain the stable value with high accuracy.
5.1.3. Identification of the Sine Load
The sampling frequency is 1000 Hz. Acceleration response of node 7 is selected for sine load identification at noise levels 1%, 3%, and 5%. The threshold is chosen as . At noise level 1%, the curves of QF and GCV and L-curve are depicted in Figures 10–12, respectively. The identified loads of the three methods are illustrated in Figure 13. The corresponding results of identified load at different noise levels are listed in Table 3.
Figures 10–13 and Table 3 show that the QFM and the L-curve method can determine the optimal regularization parameter effectively. The corresponding identified load and the real load have good correlation. Figure 11 indicates that GCV is too flat around the minimum point. In addition, the optimal regularization parameter determined by the GCV method is too large, which will result in losing too much initial information of the real response data. Thus, the accuracy of the results of sine load identification is worse. In conclusion, the adaptability of the QFM has been illustrated with different response measuring points and various types of load.
5.2. MIMO System Load Identification in Numerical Simulation
In the MIMO system situation, the same initial conditions as the SISO system situation are used. The sine load and impact load are applied on nodes 8 and 6 simultaneously. The time history of load is 0.5 s. The sine load is which is used in Section 4. A half sine wave of is used to simulate impact load. The sampling frequency is 1000 Hz. Here, the same method is used to construct the discrete linear system equation in (29). The acceleration response of nodes 5 and 7 are used to identify loads at noise levels 1%, 3%, and 5%. The threshold is selected as for the case of MIMO system. At noise level 1%, the QFM, the GCV method, and the L-curve method are used to determine the optimal regularization parameter. The curves are shown in Figures 14, 15, and 16. The corresponding identified loads are shown in Figures 17 and 18. The results of identified loads at different noise levels are illustrated in Table 4.
From Figure 14 we can see that the QFM can determine the optimal regularization parameter effectively. The results of Figures 17 and 18 and Table 4 illustrate that the identified loads which are obtained by the QFM and L-curve method have better accuracy. The optimal regularization parameter determined by the GCV method is too large, which will result in losing too much information of the real response. Thus, the accuracy of the results of GCV method is worse.
Since shape function response, which can be calculated by the finite element method, can be used to construct the system matrix. The shape function response may be polluted due to the modeling error of finite element model. It is important to verify the sensitivity of a method to a polluted system model. All other conditions being the same, modeling error is added as formula , where denotes modeling error level [35]. The results are shown in Figure 19. The QFM can determine optimal regularization parameter effectively at different modeling error levels. It can be seen that the identified results are acceptable when the modeling error level is lower than 5%. In these conditions, accuracy is not affected too much by the modeling error.
In the numerical simulation, we need several of attempts find ideal image of L-curve. However, the optimal regularization parameter can be searched by QFM directly. In addition to the impact load identification of SISO system in numerical simulation, the GCV method always has the problem that the image of GCV function is too flat and the optimal regularization parameter determined by GCV method is too large. The difference between the values of QF (quotient function) and the threshold is easier to distinguish. In addition, the QFM do not have the problem that the optimal regularization parameter is large. Therefore, the QFM can obtain the optimal regularization parameter easier and more accurately. Thus, compared with the L-curve method and GCV method, the QFM can be the best choice for determining the optimal regularization parameter.
6. The Laboratory Experiment
6.1. SISO System Load Identification in Laboratory Experiment
For the laboratory experiment of SISO system, a cantilever beam structure is considered, as shown in Figure 20. The size of the cantilever beam is 0.9 m × 0.05 m × 0.009 m. The mass density is kg/m3. The Young modulus is Gpa. The initial condition is zero. As illustrated in Figure 1, nodes are numbered 1 to 10 from left to right. Sine load is exerted at node 8 for 0.5 s whose frequency is 20 Hz. The acceleration response is measured from nodes 3 to 6 by the acceleration sensors. The sampling frequency is 1024 Hz. Here, the same method is used to construct the discrete linear system equation in (28). The threshold is chosen as . The QFM and GCV and L-curve methods are used to determine the optimal regularization parameter with the structural acceleration response of node 6 for load identification.
The curves of QF and GCV and L-curve are shown in Figures 21–23, respectively. The identification loads of these three methods are shown in Figure 24. In the laboratory experiment, the GCV function has the problems that the minimum point surrounding curve is too flat to find the minimum of the GCV function. However, the values of QF and the threshold are easy to distinguish. According to Table 5, the results validate the reasonability and effectiveness of QFM in the practice.
6.2. MIMO System Load Identification in Laboratory Experiment
For the laboratory experiment of MIMO system, the same cantilever beam structure of SISO system is considered, as shown in Figure 25. The initial condition is zero. Sine load is exerted at node 8 whose frequency is 40 Hz. Impact load is exerted at node 6 with a hammer. The acceleration response is measured at nodes 7, 9, and 10 by the acceleration sensors. The sampling frequency is 2048 Hz. Herein, the same method is used as that used in numerical simulation to construct the discrete linear system equation in (29). The threshold is chosen as . The QFM and GCV and L-curve methods are used to determine the optimal regularization parameter with the structural acceleration response of nodes 7 and 9 for load identification.
The curves of QF and GCV and L-curve are shown in Figures 26–28, respectively. The identified loads of these three methods are shown in Figures 29, 30, and 31. In the laboratory experiment of MIMO system, the GCV function also has the problem that the minimum point surrounding curve is too flat. According to the results of Table 6, the identified load corresponding to QFM has the smallest RE, which validates the reasonability and effectiveness of QFM in the practice.
From the experimental results, the L-curve method and the GCV method have the same limitations which occurred in numerical simulation. The QFM can still obtain the optimal regularization parameter easier with higher accuracy. Thus, compared with the L-curve method and the GCV method, the QFM can be the best choice for determining the optimal regularization parameter.
7. Conclusion
The numerical simulation and the laboratory experiment illustrate that it is reasonable and efficient to use QFM for determining the optimal regularization parameter in the inverse problem of load identification. The proposed method has some advantages as follows.(1)As the QFM belongs to the posterior strategy method, we do not need to find out the noise level of the structural response in advance. The method can get an approximate solution with good accuracy and stability.(2)Compared with the L-curve method, prior knowledge about the approximating range of the optimal regularization parameter is not required for the QFM. The QFM can guarantee the value of the optimal parameter within the searching range.(3)Compared with the GCV method, the difference between the values of QF and the threshold is easier to distinguish. Therefore, the optimal regularization parameter selected by the QFM can be more accurate and easier. In addition, by adjusting the value of threshold, the QFM can obtain the optimal regularization parameter which satisfies the different demands with accuracy.
Based on the above summary, the QFM can be the best choice for selecting the optimal regularization parameter. However, more sampling points will suffer more calculation burden, so finding a more efficient numerical algorithm to decrease calculated amounts is the next stage of our study.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This work is supported by the National Natural Science Foundation of China (Grant no. 11372084).