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

Authors: Wei Gao; Kaiping Yu; Ying Wu
Journal: Shock and Vibration
 (2016)

---

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


---

## Body

## 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 y
=
A
f with the discretization of Duhamel integral equation [12], where A 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], andL-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 theL-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]. TheL-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 Section2, 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 theL-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 theL-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 integral(1)
y
t
=
∫
0
t
g
t
-
τ
f
τ
d
τthe corresponding discrete linear system equation can be obtained, where y
t is a structural response which can be displacement, velocity, acceleration, and so forth. g
t is the corresponding Green kernel function which is the impulse response of the structure. f
t 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 Δ
t. We suppose there are Q time units in all and f
t
0
,
f
t
1
,
…
,
f
t
Q
-
1 are the loading points to be identified. Suppose Ω is a local support domain in the whole time domain. The local support domain Ω contains S
-
1 time units and S loading points. Without losing generality, suppose these loading points in the local support domain Ω are f
t
1
,
f
t
2
,
…
,
f
t
S and q
1
t
,
q
2
t
,
…
,
q
R
t are the basis functions. The approximate load expression can be written as(2)
f
Ω
t
=
∑
i
=
1
R
q
i
t
a
i
=
q
t
a
,where q
t
=
q
1
t
,
q
2
t
,
…
,
q
R
t is the basis function vector, and a
=
a
1
,
a
2
,
…
,
a
R
T is the combination coefficient vector. The residual between approximate load and real load in Ω is(3)
J
Ω
=
∑
j
=
1
S
f
Ω
t
j
-
f
t
j
2
=
∑
j
=
1
S
∑
i
=
1
R
q
i
t
j
a
i
-
f
t
j
2
.If we want to obtain the expression of the optimal approximate load, the number of loading points S should be equal to or larger than that of basis functions R. The necessary condition of the least residual is that the partial derivative of the combination coefficients in (3) is equal to zero. Then(4)
∂
J
∂
a
k
=
2
∑
j
=
1
S
∑
i
=
1
R
q
i
t
j
a
i
-
f
t
j
q
k
t
j
=
0
,
k
=
1,2
,
…
,
R
;then(5)
∑
j
=
1
S
∑
i
=
1
R
q
i
t
j
a
i
-
f
t
j
q
1
t
j
=
0
∑
j
=
1
S
∑
i
=
1
R
q
i
t
j
a
i
-
f
t
j
q
2
t
j
=
0
⋮
∑
j
=
1
S
∑
i
=
1
R
q
i
t
j
a
i
-
f
t
j
q
R
t
j
=
0
.Let(6)
B
=
q
1
t
1
q
1
t
2
⋯
q
1
t
S
q
2
t
1
q
2
t
2
⋯
q
2
t
S
⋮
⋮
⋱
⋮
q
R
t
1
q
R
t
2
⋯
q
R
t
S
;then(7)
a
=
a
1
,
a
2
,
…
,
a
R
T
=
B
B
T
-
1
B
f
Ω
,where f
Ω
=
f
t
1
,
f
t
2
,
…
,
f
t
S
T. Substituting (7) into (2) gives(8)
f
Ω
t
=
q
t
B
B
T
-
1
B
f
Ω
=
N
t
f
Ω
,where N
t
=
q
t
B
B
T
-
1
B 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 f
t
1 and f
t
2. The corresponding shape function vector is(9)
N
f
i
r
s
t
t
=
q
t
B
B
T
-
1
B
t
1
≤
t
≤
t
2
0
t
<
t
1
,
t
2
<
t
,where 0 denotes zero vector with the same length as the shape function vector. When the length of Δ
t is fixed, the shape function vector in (9) only depends on the basis function q
t
=
q
1
t
,
q
2
t
,
…
,
q
R
t, 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:(10)
N
~
t
=
N
~
1
t
,
N
~
2
t
,
…
,
N
~
S
t
=
N
f
i
r
s
t
t
+
t
1
=
q
t
+
t
1
B
B
T
-
1
B
0
≤
t
≤
Δ
t
0
t
<
0
,
Δ
t
<
t
.The expression of optimal approximate load in the whole time domain is(11)
f
~
h
t
=
∑
i
=
1
Q
f
~
i
h
t
=
∑
i
=
1
Q
N
~
t
-
i
-
1
Δ
t
f
~
i
,where f
~
i
h
t denotes the approximate load in i
t
h time unit and f
~
i
=
f
t
i
-
1
,
f
t
i
,
…
,
f
t
i
+
S
-
2
T. However, the loading points on the right side of f
t
Q
-
1 do not exist in practice. Without losing generality, we suppose f
t
Q
=
⋯
=
f
t
Q
+
S
-
2
=
0. Then(12)
f
~
Q
-
S
+
2
=
f
t
Q
-
S
+
1
,
…
,
f
t
Q
-
1
,
0
T
f
~
Q
-
S
+
3
=
f
t
Q
-
S
+
2
,
…
,
f
t
Q
-
1
,
0,0
T
⋮
f
~
Q
=
f
t
Q
-
1
,
0
,
…
,
0
T
.Substituting (11) into (1) gives(13)
y
t
=
∫
0
t
∑
i
=
1
Q
N
~
τ
-
i
-
1
Δ
t
f
~
i
g
t
-
τ
d
τ
.Equation (13) can be rearranged as(14)
y
t
=
∑
j
=
1
S
∑
i
=
1
Q
∫
0
t
N
~
j
τ
-
i
-
1
Δ
t
g
t
-
τ
d
τ
f
~
i
j
,where f
~
i
j denotes the j
t
h component of f
~
i. From (14), the finite element computational amount to solve the response resulting from shape function N
~
t is only S times. Suppose that, through finite element computing, the structural responses caused by each shape function load N
~
j
t are(15)
g
~
j
t
=
∑
i
=
1
Q
∫
0
t
N
~
j
τ
-
i
-
1
Δ
t
g
t
-
τ
d
τ
,
j
=
1,2
,
…
,
S
.Similar 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) gives(16)
y
=
y
t
1
,
y
t
2
,
…
,
y
t
Q
T
.With the consideration of (12) and (14), the vector of all f
~
i
j should be(17)
f
∗
=
f
t
0
,
f
t
1
,
…
,
f
t
Q
-
1
,
0
,
…
,
0
1
×
Q
+
S
-
1
T
.Similarly, the discretization of (14) gives the shape function response matrixes as(18)
G
~
j
=
0
⋯
0
0
⋯
0
⋮
⋮
⋮
0
⋯
0
︸
j
-
1
c
o
l
u
m
n
s
g
~
j
t
1
0
⋯
0
0
⋯
0
g
~
j
t
2
g
~
j
t
1
⋯
0
0
⋯
0
⋮
⋮
⋱
⋮
⋮
⋮
⋮
g
~
j
t
Q
g
~
j
t
Q
-
1
⋯
g
~
j
t
1
0
⋯
0
Q
×
Q
+
S
-
1
,
j
=
1,2
,
…
,
S
.Let(19)
G
~
=
∑
j
=
1
S
G
~
j
=
g
~
1
,
g
~
2
,
…
,
g
~
Q
+
S
-
1
,where g
~
i denotes the i
t
h column of G
~. Then we have y
=
G
~
f
∗. From (17) and (19), let(20)
G
∗
=
g
~
1
,
g
~
2
,
…
,
g
~
Q
,
(21)
f
=
f
t
0
,
f
t
1
,
…
,
f
t
Q
-
1
T
;then y
=
G
∗
f 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 supposeS
=
2
l. Then we can obtain the shape function vector corresponding to the middle and the last time unit of f
Ω
t in the local domain Ω as(22)
N
m
i
d
d
l
e
t
=
q
t
B
B
T
-
1
B
t
l
≤
t
≤
t
l
+
1
0
t
<
t
l
,
t
l
+
1
<
t
,
N
l
a
s
t
t
=
q
t
B
B
T
-
1
B
t
S
-
1
≤
t
≤
t
S
0
t
<
t
S
-
1
,
t
S
<
t
.From (9) to (11), the improved expression of optimal approximate load can be proposed as (23)
f
h
t
=
∑
i
=
1
Q
f
i
h
t
=
∑
i
=
1
l
-
1
N
~
t
-
i
-
1
Δ
t
f
~
i
+
∑
i
=
l
Q
-
l
+
1
N
-
t
-
i
-
1
Δ
t
f
-
i
+
∑
i
=
Q
-
l
+
2
Q
N
^
t
-
i
-
1
Δ
t
f
^
i
,where f
i
h
t denotes the optimal approximate load of the i
t
h time unit. Meanwhile, N
-
t
=
N
m
i
d
d
l
e
t
+
t
l, N
^
t
=
N
l
a
s
t
t
+
t
S
-
1, f
-
i
=
f
t
i
-
l
,
f
t
i
-
l
+
1
,
…
,
f
t
i
+
l
-
1
T, f
^
i
=
f
t
i
-
S
+
1
,
f
t
i
-
S
+
2
,
…
,
f
t
i
T. Without losing generality, we suppose f
t
Q
=
0. Substituting (23) into (1) gives(24)
y
t
=
∫
0
t
∑
i
=
1
l
-
1
N
~
τ
-
i
-
1
Δ
t
f
~
i
+
∑
i
=
l
Q
-
l
+
1
N
-
τ
-
i
-
1
Δ
t
f
-
i
+
∑
i
=
Q
-
l
+
2
Q
N
^
τ
-
i
-
1
Δ
t
f
^
i
g
t
-
τ
d
τ
.Similar to the construction of G
~, the shape function response matrix of N
^
t and N
-
t can be obtained with the discretization of (24) as(25)
G
-
=
g
-
1
,
g
-
2
,
…
,
g
-
Q
+
S
-
1
,
G
^
=
g
^
1
,
g
^
2
,
…
,
g
^
Q
+
S
-
1
,where g
-
i
,
g
^
i denote the i
t
h column of G
-
,
G
^. Let(26)
G
-
Q
=
g
-
11
⋯
g
-
Q
1
⋮
⋱
⋮
g
-
1
Q
-
l
+
1
⋯
g
-
Q
Q
-
l
+
1
Q
-
l
+
1
×
Q
,
G
^
Q
=
g
^
11
⋯
g
^
S
+
l
-
3
1
⋮
⋱
⋮
g
^
1
l
-
1
⋯
g
^
S
+
l
-
3
l
-
1
l
-
1
×
S
+
l
-
3
,where g
-
i
j denotes the j
t
h component of g
-
i and g
^
i
j denotes the j
t
h component of g
^
i. Let(27)
G
N
~
=
g
~
1
,
…
,
g
~
S
+
l
-
3
,
O
,
…
,
O
Q
×
Q
,
G
N
-
=
O
,
G
-
Q
T
Q
×
Q
T
,
G
N
^
=
O
O
O
G
^
Q
Q
×
Q
,where O denotes zero matrix. Let A
=
G
N
~
+
G
N
-
+
G
N
^. According to (16) and (21), the discrete linear system equation model can be obtained as(28)
y
=
A
f
.Equation (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 f
1
t and f
2
t are the two loads acting on different locations. The corresponding measured structural responses in two different locations are y
1
t and y
2
t. Therefore, according to (1), there are y
1
t
=
∫
0
t
g
11
t
-
τ
f
1
τ
+
g
21
t
-
τ
f
2
τ
d
τ and y
2
t
=
∫
0
t
g
12
t
-
τ
f
1
τ
+
g
22
t
-
τ
f
2
τ
d
τ. g
i
j
t 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 y
1
=
A
11
f
1
+
A
21
f
2 and y
2
=
A
12
f
1
+
A
22
f
2, where y
1
,
y
2
,
f
1
,
f
2 can be obtained by the discretization of y
1
t
,
y
2
t
,
f
1
t
,
f
2
t. Matrix A
i
j denotes the system matrix. Then the discrete linear system equation model of MIMO system can be constructed as(29)
y
1
y
2
=
A
11
A
21
A
21
A
22
f
1
f
2
.


## 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 as(30)
min
f
∈
R
Q
⁡
A
f
-
y
2
+
α
2
f
2
,where α
>
0 is a constant and · denotes Euler norm. The equivalent form of objective function in (30) is(31)
f
T
A
T
A
+
α
2
I
f
-
2
y
T
A
f
+
y
T
y
,where I 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 as(32)
A
T
A
+
α
2
I
f
=
A
T
y
.From SVD A
=
U
S
V
T
=
u
1
,
u
2
,
…
,
u
Q
d
i
a
g
σ
1
,
σ
2
,
…
,
σ
Q
v
1
,
v
2
,
…
,
v
Q
T and y
=
y
t
r
u
e
+
y
n
o
i
s
e [30], (32) can be rearranged as(33)
f
α
=
∑
i
=
1
Q
σ
i
2
σ
i
2
+
α
2
u
i
T
y
t
r
u
e
σ
i
v
i
+
∑
i
=
1
Q
σ
i
2
σ
i
2
+
α
2
u
i
T
y
n
o
i
s
e
σ
i
v
i
,where y
t
r
u
e is the real response data, y
n
o
i
s
e is the corresponding noise data, and σ
1
≥
σ
2
≥
⋯
≥
σ
Q
≥
0 are all the singular values of matrix A. Sometimes the system matrix A is seriously ill-conditioned. Therefore, the value corresponding to the real response y
t
r
u
e 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 y
n
o
i
s
e, which will make the error uncontrollable between the identified load and the real load. However, the regularization operator σ
i
2
/
σ
i
2
+
α
2 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 as(34)
G
C
V
α
=
y
-
A
f
α
T
r
I
-
A
A
T
A
+
α
I
-
1
A
T
,where T
r denotes the trace of matrix. The determination of optimal regularization parameter α
o
p will satisfy the condition(35)
G
C
V
α
o
p
=
min
⁡
G
C
V
α
.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]

TheL-curve criterion utilizes log-log as the scale to describe the comparison between the norm of regularization solution f
α and residual y
-
A
f
α and adopts the biggest curvature ofL-curve under the log-log scale as the criterion for determining optimal regularization parameter. Let ρ
=
l
o
g
y
-
A
f
α and θ
=
l
o
g
f
α. TheL-curve curvature is(36)
L
α
=
ρ
′
θ
′
′
-
ρ
′
′
θ
′
ρ
′
2
+
θ
′
2
3
/
2
.The condition for determining optimal regularization parameter satisfies(37)
L
α
o
p
=
max
⁡
L
α
.TheL-curve method needs to ensure the approximate range of optimal regularization parameter in advance [30]. Furthermore, there is no effective method forL-curve method to confirm optimal regularization parameter in its searching range because it is relatively small. We need several attempts to find idealL-curve [31]. However, the larger searching range will result in the overlarge calculation burden or the unsatisfactory image ofL-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 A
T
A
+
α
2
I in (31) will be well-posed. Although matrixes A and A
T
A are both ill-conditioned, the matrix A
T
A
+
α
2
I is well-posed when the value of regularization parameter α is large enough. As a result, a good least square solution f
α that satisfies (32) can be obtained as(38)
A
T
y
-
A
T
A
f
α
=
α
2
f
α
;then (38) can be rearranged as(39)
A
T
y
-
A
T
A
f
α
=
α
2
f
α
.From (33), f
α 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 as(40)
A
T
y
-
A
T
A
f
α
α
2
f
α
=
1
.With the value of regularization parameter α decreasing, matrix A
T
A
+
α
2
I will transform from being well-posed to being weakly ill-conditioned. Thus the least square solution f
α that is solved through the quadratic programming algorithm is not fairly accurate. Then (40) becomes(41)
A
T
y
-
A
T
A
f
α
α
2
f
α
≈
1
.If the value of regularization parameter α decreases more, matrix A
T
A
+
α
2
I will transform to being strongly ill-conditioned. The accuracy of the least square solution f
α will become even worse. Then (41) becomes(42)
A
T
y
-
A
T
A
f
α
α
2
f
α
≠
1
.Therefore, quotient function is defined as(43)
H
α
=
A
T
y
-
A
T
A
f
α
α
2
f
α
.We hope to find out a parameter α
0
>
0 satisfying the following conditions:(i)
ifα
≥
α
0, then H
α
=
A
T
y
-
A
T
A
f
α
/
α
2
f
α
=
1;
(ii)
ifα
<
α
0, then H
α
=
A
T
y
-
A
T
A
f
α
/
α
2
f
α
≠
1.Then the parameter α
0
>
0 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σ
1 by performing SVD to the system matrix A. The initial reference parameter vector can be defined as β
1
=
σ
1
/
100
,
2
σ
1
/
100
,
…
,
σ
1.
(2)
Calculate values of the QF corresponding to every component in the vectorβ
1. The vector that is made of the values of QF can be obtained as H
=
H
σ
1
/
100
,
H
2
σ
1
/
100
,
…
,
H
σ
1.
(3)
Set a thresholde
o
p that equals 1 approximately. All values of QF in vector H are checked. If the distance between any component and 1 is less than or equal to e
o
p, then the reference parameter vector can be renewed as β
2
=
β
1
/
2.
(4)
Repeat steps (2) and (3) until the distance between a certain component of vector H and 1 is more than e
o
p (suppose the component is the j
-
1
t
h component of H). Now we suppose the reference parameter vector is β
i
=
β
1
/
2
i
-
1
=
β
i
1
,
β
i
2
,
…
,
β
i
100.
(5)
Stop the searching process with the optimal regularization parameter being thej
t
h component of β
i; that is, α
o
p
=
β
i
j.The algorithm realization of the QFM can be explained as follows: the optimal regularization parameterα
o
p is smaller than the largest singular value σ
1. Thus, the largest component of the initial reference vector β
1 is chosen as the largest singular value σ
1, which can guarantee the value of the optimal regularization parameter α
o
p within the searching range. According to the different demands for accuracy in practice application, the threshold e
o
p can be chosen as 10
-
3, 10
-
4, or 10
-
5. If there is no demand, we choose threshold e
o
p as 10
-
4. According to the threshold e
o
p, the searching process will be stopped if the distance between any component of vector H and 1 is more than e
o
p. Then we can determine the optimal regularization parameter α
o
p from the conditions (i) and (ii) in Section 4.3.


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

According to the discrete linear system equation in (28), the GCV function is defined as(34)
G
C
V
α
=
y
-
A
f
α
T
r
I
-
A
A
T
A
+
α
I
-
1
A
T
,where T
r denotes the trace of matrix. The determination of optimal regularization parameter α
o
p will satisfy the condition(35)
G
C
V
α
o
p
=
min
⁡
G
C
V
α
.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]

TheL-curve criterion utilizes log-log as the scale to describe the comparison between the norm of regularization solution f
α and residual y
-
A
f
α and adopts the biggest curvature ofL-curve under the log-log scale as the criterion for determining optimal regularization parameter. Let ρ
=
l
o
g
y
-
A
f
α and θ
=
l
o
g
f
α. TheL-curve curvature is(36)
L
α
=
ρ
′
θ
′
′
-
ρ
′
′
θ
′
ρ
′
2
+
θ
′
2
3
/
2
.The condition for determining optimal regularization parameter satisfies(37)
L
α
o
p
=
max
⁡
L
α
.TheL-curve method needs to ensure the approximate range of optimal regularization parameter in advance [30]. Furthermore, there is no effective method forL-curve method to confirm optimal regularization parameter in its searching range because it is relatively small. We need several attempts to find idealL-curve [31]. However, the larger searching range will result in the overlarge calculation burden or the unsatisfactory image ofL-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 A
T
A
+
α
2
I in (31) will be well-posed. Although matrixes A and A
T
A are both ill-conditioned, the matrix A
T
A
+
α
2
I is well-posed when the value of regularization parameter α is large enough. As a result, a good least square solution f
α that satisfies (32) can be obtained as(38)
A
T
y
-
A
T
A
f
α
=
α
2
f
α
;then (38) can be rearranged as(39)
A
T
y
-
A
T
A
f
α
=
α
2
f
α
.From (33), f
α 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 as(40)
A
T
y
-
A
T
A
f
α
α
2
f
α
=
1
.With the value of regularization parameter α decreasing, matrix A
T
A
+
α
2
I will transform from being well-posed to being weakly ill-conditioned. Thus the least square solution f
α that is solved through the quadratic programming algorithm is not fairly accurate. Then (40) becomes(41)
A
T
y
-
A
T
A
f
α
α
2
f
α
≈
1
.If the value of regularization parameter α decreases more, matrix A
T
A
+
α
2
I will transform to being strongly ill-conditioned. The accuracy of the least square solution f
α will become even worse. Then (41) becomes(42)
A
T
y
-
A
T
A
f
α
α
2
f
α
≠
1
.Therefore, quotient function is defined as(43)
H
α
=
A
T
y
-
A
T
A
f
α
α
2
f
α
.We hope to find out a parameter α
0
>
0 satisfying the following conditions:(i)
ifα
≥
α
0, then H
α
=
A
T
y
-
A
T
A
f
α
/
α
2
f
α
=
1;
(ii)
ifα
<
α
0, then H
α
=
A
T
y
-
A
T
A
f
α
/
α
2
f
α
≠
1.Then the parameter α
0
>
0 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σ
1 by performing SVD to the system matrix A. The initial reference parameter vector can be defined as β
1
=
σ
1
/
100
,
2
σ
1
/
100
,
…
,
σ
1.
(2)
Calculate values of the QF corresponding to every component in the vectorβ
1. The vector that is made of the values of QF can be obtained as H
=
H
σ
1
/
100
,
H
2
σ
1
/
100
,
…
,
H
σ
1.
(3)
Set a thresholde
o
p that equals 1 approximately. All values of QF in vector H are checked. If the distance between any component and 1 is less than or equal to e
o
p, then the reference parameter vector can be renewed as β
2
=
β
1
/
2.
(4)
Repeat steps (2) and (3) until the distance between a certain component of vector H and 1 is more than e
o
p (suppose the component is the j
-
1
t
h component of H). Now we suppose the reference parameter vector is β
i
=
β
1
/
2
i
-
1
=
β
i
1
,
β
i
2
,
…
,
β
i
100.
(5)
Stop the searching process with the optimal regularization parameter being thej
t
h component of β
i; that is, α
o
p
=
β
i
j.The algorithm realization of the QFM can be explained as follows: the optimal regularization parameterα
o
p is smaller than the largest singular value σ
1. Thus, the largest component of the initial reference vector β
1 is chosen as the largest singular value σ
1, which can guarantee the value of the optimal regularization parameter α
o
p within the searching range. According to the different demands for accuracy in practice application, the threshold e
o
p can be chosen as 10
-
3, 10
-
4, or 10
-
5. If there is no demand, we choose threshold e
o
p as 10
-
4. According to the threshold e
o
p, the searching process will be stopped if the distance between any component of vector H and 1 is more than e
o
p. Then we can determine the optimal regularization parameter α
o
p 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 Section6, as shown in Figure 20. From left to right are nodes 1 to 10, as shown in Figure 1. Take f
1
=
20
sin
⁡
60
π
t
1
-
c
o
s
20
π
t, f
2
=
40
sin
⁡
40
π
t, and f
3
=
40
sin
⁡
200
π
t as the external load applied on node 8. Amplitude modulation (AM) load f
1 and sine load f
2 are acting on node 8 for 0.3 s. Sine load f
3 is acting on node 8 for 0.005 s to simulate impact load.Figure 1
Node location of the cantilever beam model.In the numerical examples,q
t
=
[
1
,
t
,
t
2
] is selected as the basis function vector to construct the system matrix in (28). The number of loading points S
=
4. We can solve (28) by using the Tikhonov regularization method and determine the optimal regularization parameter through the QFM and the GCV and theL-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 as(44)
y
=
y
t
r
u
e
+
l
n
o
n
o
i
s
e
s
t
d
y
t
r
u
e
,where y is the acceleration response for load identification; y
t
r
u
e is the numerical calculated acceleration response; l
n
o is the noise level; n
o
i
s
e is a vector of white noise with zero mean and unit standard deviation; s
t
d
· denotes the standard deviation. The relative error (RE) and correlation coefficient (CC) between the identified load and the real load can be defined as(45)
R
E
=
f
α
-
f
r
e
a
l
f
r
e
a
l
×
100
%
,
C
C
=
∑
k
=
1
Q
f
α
k
-
E
f
α
f
r
e
a
l
k
-
E
f
r
e
a
l
f
α
-
E
f
α
f
r
e
a
l
-
E
f
r
e
a
l
,where f
α
k denotes the k
t
h component of f
α, f
r
e
a
l denotes the real load, and E
· 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 ase
o
p
=
1
0
-
4. 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 α
<
0.0015, the values of QF are far away from 1. When α
≥
0.0015, 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 0.0015.Figure 2
QF for AM load identification of SISO system.Similarly, the regularization parameter is0.0024 based on the GCV method as shown in Figure 3. The result ofL-curve method is 1.3310
×
1
0
-
5 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.Table 1
Relative results of identified AM load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0015
0.0024
1.331× 10−5
3.34
3.34
2.83
0.9996
0.9932
0.9997
3
0.0020
0.0053
1.781× 10−5
7.63
7.54
6.84
0.9973
0.9974
0.9979
5
0.0025
0.0069
1.391× 10−5
12.17
11.98
11.13
0.9930
0.9932
0.9942Figure 3
GCV for AM load identification of SISO system.Figure 4
L-curve for AM load identification of SISO system.Figure 5
Identified AM load of SISO system.


#### 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 ase
o
p
=
1
0
-
4. 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 andL-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.Table 2
Relative results of identified impact load of SISO system.
Real (N)
Noise level (%)
1
3
5
Method
IIL (N)
RE (%)
IIL (N)
RE (%)
IIL (N)
RE (%)
Amplitude = 39.92
QFM
40.03
0.26
41.36
3.62
40.69
1.94
GCV
38.67
3.14
39.39
1.33
39.61
0.76
L-curve
40.19
0.68
41.62
4.25
43.04
7.81Figure 6
QF for impact load identification of SISO system.Figure 7
GCV for impact load identification of SISO system.Figure 8
L-curve for impact load identification of SISO system.Figure 9
Identified impact load of SISO system.Figures6 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 ase
o
p
=
1
0
-
4. At noise level 1%, the curves of QF and GCV andL-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.Table 3
Relative results of identified sine load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0026
0.2110
0.0032
4.02
12.55
4.0
0.9996
0.9999
0.9996
3
0.0026
0.2110
0.0028
9.32
12.87
9.32
0.9969
0.9992
0.9969
5
0.0621
0.2110
0.0023
12.54
13.58
14.91
0.9927
0.9979
0.9913Figure 10
QF for sine load identification of SISO system.Figure 11
GCV for sine load identification of SISO system.Figure 12
L-curve for sine load identification of SISO system.Figure 13
Identified sine load of SISO system.Figures10–13 and Table 3 show that the QFM and theL-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 isf
2 which is used in Section 4. A half sine wave of f
4
t
=
45
sin
⁡
100
π
t 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 e
o
p
=
1
0
-
4 for the case of MIMO system. At noise level 1%, the QFM, the GCV method, and theL-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.Table 4
Relative results of identified load of MIMO system.
Noise level (%)
QFM
GCV
L-curve
α
op
RE (%)
α
op
RE (%)
α
op
RE (%)
f
2
f
4
f
2
f
4
f
2
f
4
1
0.0031
7.45
6.80
0.1900
22.51
20.75
0.0032
7.45
6.80
3
0.0047
8.84
6.34
0.1900
22.64
20.00
0.0037
8.82
6.28
5
0.0189
11.56
7.21
0.1800
22.48
18.57
0.0041
11.07
5.75Figure 14
QF for load identification of MIMO system.Figure 15
GCV for load identification of MIMO system.Figure 16
L-curve for load identification of MIMO system.Figure 17
Identified sine load of MIMO system.Figure 18
Identified impact load of MIMO system.From Figure14 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 andL-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 formulaA
e
r
r
o
r
=
A
1
+
Δ, 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.Figure 19
Results at different modeling error levels.Figure 20
The cantilever beam structure of SISO system.In the numerical simulation, we need several of attempts find ideal image ofL-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 theL-curve method and GCV method, the QFM can be the best choice for determining the optimal regularization parameter.


## 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 ase
o
p
=
1
0
-
4. 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 α
<
0.0015, the values of QF are far away from 1. When α
≥
0.0015, 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 0.0015.Figure 2
QF for AM load identification of SISO system.Similarly, the regularization parameter is0.0024 based on the GCV method as shown in Figure 3. The result ofL-curve method is 1.3310
×
1
0
-
5 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.Table 1
Relative results of identified AM load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0015
0.0024
1.331× 10−5
3.34
3.34
2.83
0.9996
0.9932
0.9997
3
0.0020
0.0053
1.781× 10−5
7.63
7.54
6.84
0.9973
0.9974
0.9979
5
0.0025
0.0069
1.391× 10−5
12.17
11.98
11.13
0.9930
0.9932
0.9942Figure 3
GCV for AM load identification of SISO system.Figure 4
L-curve for AM load identification of SISO system.Figure 5
Identified AM load of SISO system.


### 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 ase
o
p
=
1
0
-
4. 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 andL-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.Table 2
Relative results of identified impact load of SISO system.
Real (N)
Noise level (%)
1
3
5
Method
IIL (N)
RE (%)
IIL (N)
RE (%)
IIL (N)
RE (%)
Amplitude = 39.92
QFM
40.03
0.26
41.36
3.62
40.69
1.94
GCV
38.67
3.14
39.39
1.33
39.61
0.76
L-curve
40.19
0.68
41.62
4.25
43.04
7.81Figure 6
QF for impact load identification of SISO system.Figure 7
GCV for impact load identification of SISO system.Figure 8
L-curve for impact load identification of SISO system.Figure 9
Identified impact load of SISO system.Figures6 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 ase
o
p
=
1
0
-
4. At noise level 1%, the curves of QF and GCV andL-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.Table 3
Relative results of identified sine load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0026
0.2110
0.0032
4.02
12.55
4.0
0.9996
0.9999
0.9996
3
0.0026
0.2110
0.0028
9.32
12.87
9.32
0.9969
0.9992
0.9969
5
0.0621
0.2110
0.0023
12.54
13.58
14.91
0.9927
0.9979
0.9913Figure 10
QF for sine load identification of SISO system.Figure 11
GCV for sine load identification of SISO system.Figure 12
L-curve for sine load identification of SISO system.Figure 13
Identified sine load of SISO system.Figures10–13 and Table 3 show that the QFM and theL-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.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 ase
o
p
=
1
0
-
4. 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 α
<
0.0015, the values of QF are far away from 1. When α
≥
0.0015, 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 0.0015.Figure 2
QF for AM load identification of SISO system.Similarly, the regularization parameter is0.0024 based on the GCV method as shown in Figure 3. The result ofL-curve method is 1.3310
×
1
0
-
5 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.Table 1
Relative results of identified AM load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0015
0.0024
1.331× 10−5
3.34
3.34
2.83
0.9996
0.9932
0.9997
3
0.0020
0.0053
1.781× 10−5
7.63
7.54
6.84
0.9973
0.9974
0.9979
5
0.0025
0.0069
1.391× 10−5
12.17
11.98
11.13
0.9930
0.9932
0.9942Figure 3
GCV for AM load identification of SISO system.Figure 4
L-curve for AM load identification of SISO system.Figure 5
Identified AM load of SISO system.


## 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 ase
o
p
=
1
0
-
4. 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 andL-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.Table 2
Relative results of identified impact load of SISO system.
Real (N)
Noise level (%)
1
3
5
Method
IIL (N)
RE (%)
IIL (N)
RE (%)
IIL (N)
RE (%)
Amplitude = 39.92
QFM
40.03
0.26
41.36
3.62
40.69
1.94
GCV
38.67
3.14
39.39
1.33
39.61
0.76
L-curve
40.19
0.68
41.62
4.25
43.04
7.81Figure 6
QF for impact load identification of SISO system.Figure 7
GCV for impact load identification of SISO system.Figure 8
L-curve for impact load identification of SISO system.Figure 9
Identified impact load of SISO system.Figures6 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 ase
o
p
=
1
0
-
4. At noise level 1%, the curves of QF and GCV andL-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.Table 3
Relative results of identified sine load of SISO system.
Noise level (%)
α
op
RE (%)
CC
QFM
GCV
L-curve
QFM
GCV
L-curve
QFM
GCV
L-curve
1
0.0026
0.2110
0.0032
4.02
12.55
4.0
0.9996
0.9999
0.9996
3
0.0026
0.2110
0.0028
9.32
12.87
9.32
0.9969
0.9992
0.9969
5
0.0621
0.2110
0.0023
12.54
13.58
14.91
0.9927
0.9979
0.9913Figure 10
QF for sine load identification of SISO system.Figure 11
GCV for sine load identification of SISO system.Figure 12
L-curve for sine load identification of SISO system.Figure 13
Identified sine load of SISO system.Figures10–13 and Table 3 show that the QFM and theL-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 isf
2 which is used in Section 4. A half sine wave of f
4
t
=
45
sin
⁡
100
π
t 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 e
o
p
=
1
0
-
4 for the case of MIMO system. At noise level 1%, the QFM, the GCV method, and theL-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.Table 4
Relative results of identified load of MIMO system.
Noise level (%)
QFM
GCV
L-curve
α
op
RE (%)
α
op
RE (%)
α
op
RE (%)
f
2
f
4
f
2
f
4
f
2
f
4
1
0.0031
7.45
6.80
0.1900
22.51
20.75
0.0032
7.45
6.80
3
0.0047
8.84
6.34
0.1900
22.64
20.00
0.0037
8.82
6.28
5
0.0189
11.56
7.21
0.1800
22.48
18.57
0.0041
11.07
5.75Figure 14
QF for load identification of MIMO system.Figure 15
GCV for load identification of MIMO system.Figure 16
L-curve for load identification of MIMO system.Figure 17
Identified sine load of MIMO system.Figure 18
Identified impact load of MIMO system.From Figure14 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 andL-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 formulaA
e
r
r
o
r
=
A
1
+
Δ, 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.Figure 19
Results at different modeling error levels.Figure 20
The cantilever beam structure of SISO system.In the numerical simulation, we need several of attempts find ideal image ofL-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 theL-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 Figure20. The size of the cantilever beam is 0.9 m × 0.05 m × 0.009 m. The mass density is ρ
=
7.8
×
10
3 kg/m3. The Young modulus is E
=
200 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 e
o
p
=
1
0
-
4. The QFM and GCV andL-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 andL-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.Table 5
Relative results of laboratory experiment of SISO system.
Method
α
op
RE (%)
CC
QFM
0.0642
8.75
0.9983
GCV
0.1270
11.62
0.9983
L-curve
1.9010× 10−7
8.01
0.9982Figure 21
QF in experiment of SISO system.Figure 22
GCV in experiment of SISO system.Figure 23
L-curve in experiment of SISO system.Figure 24
Sine loads in experiment of SISO system.


### 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 Figure25. 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 e
o
p
=
1
0
-
4. The QFM and GCV andL-curve methods are used to determine the optimal regularization parameter with the structural acceleration response of nodes 7 and 9 for load identification.Figure 25
The cantilever beam structure of MIMO system.The curves of QF and GCV andL-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.Table 6
Relative results of laboratory experiment of MIMO system.
Method
Sine load
Impact load
α
op
RE (%)
Real (N)
IIL (N)
RE (%)
QFM
0.0337
4.63
33.0
31.31
5.03
GCV
0.1880
10.09
25.38
23.10
L-curve
9.00× 10−6
5.28
31.34
5.05Figure 26
QF in experiment of MIMO system.Figure 27
GCV in experiment of MIMO system.Figure 28
L-curve in experiment of MIMO system.Figure 29
Sine load in experiment of MIMO system.Figure 30
Impact load in experiment of MIMO system.Figure 31
Local impact load in experiment of MIMO system.From the experimental results, theL-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 theL-curve method and the GCV method, the QFM can be the best choice for determining the optimal regularization parameter.


## 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 Figure20. The size of the cantilever beam is 0.9 m × 0.05 m × 0.009 m. The mass density is ρ
=
7.8
×
10
3 kg/m3. The Young modulus is E
=
200 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 e
o
p
=
1
0
-
4. The QFM and GCV andL-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 andL-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.Table 5
Relative results of laboratory experiment of SISO system.
Method
α
op
RE (%)
CC
QFM
0.0642
8.75
0.9983
GCV
0.1270
11.62
0.9983
L-curve
1.9010× 10−7
8.01
0.9982Figure 21
QF in experiment of SISO system.Figure 22
GCV in experiment of SISO system.Figure 23
L-curve in experiment of SISO system.Figure 24
Sine loads in experiment of SISO system.


## 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 Figure25. 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 e
o
p
=
1
0
-
4. The QFM and GCV andL-curve methods are used to determine the optimal regularization parameter with the structural acceleration response of nodes 7 and 9 for load identification.Figure 25
The cantilever beam structure of MIMO system.The curves of QF and GCV andL-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.Table 6
Relative results of laboratory experiment of MIMO system.
Method
Sine load
Impact load
α
op
RE (%)
Real (N)
IIL (N)
RE (%)
QFM
0.0337
4.63
33.0
31.31
5.03
GCV
0.1880
10.09
25.38
23.10
L-curve
9.00× 10−6
5.28
31.34
5.05Figure 26
QF in experiment of MIMO system.Figure 27
GCV in experiment of MIMO system.Figure 28
L-curve in experiment of MIMO system.Figure 29
Sine load in experiment of MIMO system.Figure 30
Impact load in experiment of MIMO system.Figure 31
Local impact load in experiment of MIMO system.From the experimental results, theL-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 theL-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 theL-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.



---

*Source: 7328969-2016-03-16.xml*