Accelerated generalized minimax-concave sparse regularization for impact force reconstruction and localization

Abstract

Impact force identification has always been of significance for structure health monitoring especially on the applications involving composite materials. As a typical inverse problem, impact force reconstruction and localization is undoubtedly a challenging task. The well-known ℓ₁ sparse regularization has a tendency to underestimate the amplitude of impact forces. To alleviate this limitation, we propose an accelerated generalized minimax-concave (AGMC) for sparse regularization that employs a non-convex generalized minimax-concave (GMC) penalty as the regularizer and incorporates an acceleration technique to expedite the attainment of the global minimum. Compared with the classic ℓ₁-norm penalty, the GMC penalty can not only induce sparsity in the estimation, but also maintain the convexity of the cost function, so that the global optimal solution can be obtained through convex optimization algorithms. This method is applied to solve the impact force identification problem with unknown force locations to simultaneously reconstruct and localize impact forces in the under-determined case utilizing a limited number of sensors. Meanwhile, K-sparsity criterion is used to adaptively select regularization parameters by taking advantage of the sparse prior knowledge on impact forces. Simulations and experiments are conducted on a composite plate to verify the computational efficiency and robustness of the AGMC method in terms of impact force reconstruction and localization, particularly in the presence of noise. Results demonstrate that the proposed AGMC method achieves faster convergence and provides more accurate and sparse reconstruction and localization of impact forces compared to other state-of-the-art sparse regularization methods.

Introduction

Composite materials are widely used in the mechanical engineering fields, especially in aeronautical structures due to their excellent properties such as high specific stiffness and specific strength. ¹ However, they have the disadvantage of poor impact resistance. Composite structures are susceptible to barely visible impact damage (BVID) such as debonding and delamination when impacted by foreign objects such as birds, rocks, and hail. ² If left undetected, BVIDs can pose a significant threat to the integrity of the structure, potentially leading to operational issues and irreversible damage. ³ Impact force identification is an important technique that can aid in quickly identifying where impact damage has occurred and help assess the structural integrity of composite structures. Therefore, identifying impact forces acting on composite structures, including reconstructing their time histories and localizing impact positions, is crucial for assessing the structural integrity and detecting any potential damage. Considering the impracticality of directly measuring impact forces at unknown locations using force sensors, measurable structural responses are instead utilized for inverse solving of the impact forces. Then, many methods have been developed for addressing this typical inverse problem of impact force identification. ⁴

To mitigate the highly under-determined and unstable nature of the highly ill-conditioned inverse problem of impact force identification, regularization techniques are commonly utilized. ⁵ Tikhonov regularization method is one of the classical methods for the inverse problem. Jacquelin et al. ⁶ compared generalized singular value decomposition (GSVD), Tikhonov, and TSVD methods to reconstruct impact forces imposed on an aluminum plate in time domain, and stated that the condition number of the transfer function is influenced by the position of measuring points. Li et al. ⁷ employed the Nelder–Mead method for force localization and subsequently utilized Tikhonov regularization for the reconstruction of their time histories in a cantilever beam with two accelerometers. Yan et al. ⁸ realized impact force identification with two nested loops, that is, an outer loop to localize the impact force with a non-linear unscented Kalman filter, and an inner loop to reconstruct the time-history with Tikhonov regularization. Although ℓ₂ regularization methods are straightforward and practical, their performance for impact force identification is typically unsatisfactory when dealing with under-determined sensor placement cases.

As the research interest in sparse regularization methods grows, some scholars have shifted their attention from traditional ℓ₂-norm-based methods to solving cost functions utilizing the ℓ₁ norm. This is due to the fact that ℓ₁ regularization methods encourage solutions with a minimum number of non-zero values, which aligns with the sparse nature of impact forces in the joint time-space domain, thus making these methods to work well in the case of under-determined sensor configurations. ⁹ Since the ℓ₁-norm minimization model is convex, it can be solved by convex algorithms, such as the gradient projection method, ¹⁰ the interior point method, ¹¹ the stage-wise orthogonal matching pursuit method, ¹² and the iterative soft threshold algorithm (ISTA). ¹³ Ginsberg et al. ¹⁴ proposed a method for simultaneous identification of impact locations and time histories of impact forces on a simple beam structure using direct deconvolution, which involved solving an extended ℓ₁-minimization problem. Rezayat et al. ¹⁵ identified point forces imposed on a footbridge via a regularization method that combines ℓ₁ and ℓ₂-norms. Pan et al. ¹⁶ proposed a weighted ℓ₁-norm regularization, and it performed better than the ℓ₁-norm and Tikhonov regularization methods in identifying the moving force acting on a steel beam. Aucejo et al. ¹⁷ proposed a space-frequency multiplicative regularization method considering the sparsity of excitation sources in space domain and verified this method on a thin simply supported steel beam. Qiao et al. ¹⁸ proposed an enhanced sparse regularization method on the basis of reweighted ℓ₁-norm minimization for impact force identification. While ℓ₁ norm has been widely adopted for promoting sparsity in the solution of impact force identification problems, its effectiveness in achieving accurate estimates of impact forces is often limited and may lead to underestimation of force amplitudes. ¹⁹

In recent years, researchers have investigated the use of non-convex penalties in the impact force identification problem as a means of addressing the limitations of ℓ₁ regularization and promoting sparsity. Chartrand et al. ²⁰ indicated that non-convex functions such as ℓ _p quasi-norm (0 < p < 1) exhibit superior sparsity-inducing capabilities compared to convex ℓ₁ regularization, especially in cases where the number of measurements is limited. Aucejo et al. ²¹ proposed the incorporation of a local regularization term, $\mathcal{R}\left({\mathbf{F}}_{\mathbf{r}}\right)={\Vert {\mathbf{F}}_{\mathbf{r}}\Vert }_{q_r}^{q_r}$ , into the minimization problem, where the sparsity of $\mathcal{R}\left({\mathbf{F}}_{\mathbf{r}}\right)$ is determined by the value of q _r which is set to less than 1. Liu et al. ²² developed a non-convex penalty-based optimization objective function for impact force identification, which was solved using a fast non-convex overlapping group sparsity algorithm and validated on a stiffened composite plate. After Selesnick ²³ proposed a non-convex penalty named generalized minimax-concave (GMC) for de-noising without underestimating high-amplitude components of analyzed signals, this non-convex penalty has gradually came into view. As Ref. ²³ emphasized, through introducing the non-convex GMC penalty into the cost function, the convexity of the cost function can still be guaranteed under a certain condition, so that the function can be solved by convex optimization algorithms. From then on, the GMC penalty was widely used for signal processing in bearing fault diagnosis,^24,25 gearbox fault diagnosis,^26,27 ship wake detection, ²⁸ acoustic source reconstruction, ²⁹ and so forth. In the field of impact force identification, Liu et al. ³⁰ introduced the GMC penalty into the objective function and verified in simulations and experiments that the GMC method outperforms ℓ₁ regularization in reconstructing time histories of impact forces with known force locations using one accelerometer under the even-determined condition.

It is worth noting that in practical engineering problems, the locations of impact forces cannot be known in advance. Thus, the primary objective of this contribution is to address the challenge of monitoring a greater number of impact locations with a reduced number of sensors. This requires solving a large-scale under-determined inverse problem, which is a significant research problem in the field. Although the GMC method can considerably improve the accuracy of the amplitude estimation, it will suffer from a slow convergence when dealing with such a large-scale problem with a high-dimension transfer matrix. To the best of authors’ knowledge, the GMC method has only been used to solve the impact force identification problem in the even-determined case so far, as in Ref. ³⁰ the application of the GMC method in the under-determined case has not been investigated. Therefore, we propose a novel impact force identification method named accelerated generalized minimax-concave sparse regularization (AGMC) to simultaneously reconstruct and localize impact forces at unknown impact locations in the under-determined case, which can provide a faster convergence rate and higher identification accuracy under the same iteration termination criteria. This method utilizes an acceleration strategy in solving the cost function including the non-convex GMC penalty, which not only inherits the advantages of the GMC method but also accelerates the convergence speed and improves the computational efficiency. The K-sparsity criterion is used to select the regularization parameter of this method.

The rest of this paper is organized as follows. The Dynamic modeling section briefly deduces the modeling problem for impact force identification. The Sparse regularization methods for impact force reconstruction and localization section describes the existing ISTA and GMC methods, derives the AGMC method, and uses the K-sparsity criterion to uniformly select regularization parameters for the three regularization methods. The Numerical verification section carries out numerical verification on a composite plate. The Experimental verification section verifies the AGMC method in an experiment on a composite plate with the similar parameters and constraints as the simulation. The Conclusions section concludes this paper.

Dynamic modeling

In the context of a linear time-invariant single-input single-output dynamic system, the convolution operation can be used to express the relationship between the input excitation and output response, which is

$s\left(t\right)=a\left(t\right)\otimes f\left(t\right)={\int }_0^{t}a\left(t-\tau \right)f\left(t\right)d\tau$

(1)

where s(t) denotes the system response, a(t) represents the impulse response function, f(t) is the excitation impact force, and the symbol ⊗ represents the convolution operation. It is assumed that s(t) = a(t) = f(t) = 0 when t < 0.

The continuous convolution model as presented in equation (1) can be discretized as

$s\left(N\mathrm{\Delta }t\right)=\mathrm{\Delta }t{\displaystyle \sum }_{i=1}^{N}a\left(\left(N-i\right)\mathrm{\Delta }t\right)f\left(i\mathrm{\Delta }t\right)$

(2)

where N denotes the data length of the discretized impulse response function(IRF), and Δt represents the sampling interval. Equation (2) can be expressed in an alternative form as

$\left[\begin{array}{c}s\left(\mathrm{\Delta }t\right)\\ s\left(2\mathrm{\Delta }t\right)\\ ⋮\\ s\left(\left(N-1\right)\mathrm{\Delta }t\right)\\ s\left(N\mathrm{\Delta }t\right)\end{array}\right]=\mathrm{\Delta }t\left[\begin{array}{ccccc}a\left(\mathrm{\Delta }t\right)& 0& \cdots & 0& 0\\ a\left(2\mathrm{\Delta }t\right)& a\left(\mathrm{\Delta }t\right)& \cdots & 0& 0\\ ⋮& ⋮& \cdots & ⋮& ⋮\\ a\left(\left(N-1\right)\mathrm{\Delta }t\right)& a\left(\left(N-2\right)\mathrm{\Delta }t\right)& \cdots & a\left(\mathrm{\Delta }t\right)& 0\\ a\left(N\mathrm{\Delta }t\right)& a\left(\left(N-1\right)\mathrm{\Delta }t\right)& \cdots & a\left(2\mathrm{\Delta }t\right)& a\left(\mathrm{\Delta }t\right)\end{array}\right]\left[\begin{array}{c}f\left(\mathrm{\Delta }t\right)\\ f\left(2\mathrm{\Delta }t\right)\\ ⋮\\ f\left(\left(N-1\right)\mathrm{\Delta }t\right)\\ f\left(N\mathrm{\Delta }t\right)\end{array}\right]$

(3)

For the purpose of brevity and convenience, equation (3) can be compactly represented as

${\mathbf{s}}_{\mathrm{i}}={\mathbf{A}}_{\text{ij}}{\mathbf{f}}_{\mathrm{j}}$

(4)

where the vector ${\mathbf{s}}_{\mathrm{i}}\in {\mathfrak{R}}^N$ is the response at a certain position i, the vector ${\mathbf{f}}_{\mathrm{j}}\in {\mathfrak{R}}^N$ is the excitation force at a certain position j, and ${\mathbf{A}}_{\text{ij}}\in {\mathfrak{R}}^{N\times N}$ refers to the transfer matrix between the excitation position j and the response position i.

When impact forces are applied to various unknown locations on a structure, the corresponding responses from multiple positions are simultaneously recorded as

$\left[\begin{array}{l}{\mathbf{s}}_1\\ {\mathbf{s}}_2\\ ⋮\\ {\mathbf{s}}_m\end{array}\right]=\left[\begin{array}{llll}{\mathbf{A}}_{11}& {\mathbf{A}}_{12}& \cdots & {\mathbf{A}}_{1n}\\ {\mathbf{A}}_{21}& {\mathbf{A}}_{22}& \cdots & {\mathbf{A}}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{A}}_{m1}& {\mathbf{A}}_{m2}& \cdots & {\mathbf{A}}_{mn}\end{array}\right]\left[\begin{array}{l}{\mathbf{f}}_1\\ {\mathbf{f}}_2\\ ⋮\\ {\mathbf{f}}_n\end{array}\right]$

(5)

where m and n denotes the respective quantities of measurement responses and impact force excitations.

Equation (5) can also be concisely expressed as the following compact form

$\mathbf{s}=\mathbf{A}\mathbf{f}+\mathbf{e}$

(6)

where the response vector $\mathbf{s}\in {\mathfrak{R}}^{mN}$ , the force vector $\mathbf{f}\in {\mathfrak{R}}^{nN}$ , the block Toeplitz-like matrix $\mathbf{A}\in {\mathfrak{R}}^{mN\times nN}$ is commonly referred to as the transfer matrix in literature, and the vector $\mathbf{e}\in {\mathfrak{R}}^{mN}$ corresponds to the measurement noise that is inherent in the actual measurement process. Based on the relationship between the number of measurements m and excitations n, the inverse problem of impact force identification can be classified into three distinct categories: (I) if m > n, the inverse problem is in an over-determined condition; (II) if m = n, the inverse problem is in an even-determined condition; (III) if m < n, the inverse problem is in an under-determined condition. Given that it is often unrealistic for the number of sensors m to exceed the number of potential impact locations n in real-world scenarios, ³¹ this paper focuses on investigating the under-determined inverse problem. By solving this problem, it is possible to achieve both the localization and time-history reconstruction of impact forces.

Sparse regularization methods for impact force reconstruction and localization

The impact force identification problem is a typical ill-posed problem. The condition number of the transfer matrix A is quite large, so the inevitable noise in the response will cause a large error in the solution result. Therefore, regularization methods are usually applied to solve this problem. Since the traditional ℓ₂ regularization methods fail to solve this problem in the under-determined case, ³² sparse regularization methods are developed using the sparsity of impact forces in time-space domain as prior knowledge to improve the accuracy of identification results. In this section, we will offer a brief overview of the ℓ₁ sparse regularization method and the generalized minimax-concave (GMC) sparse regularization method. Furthermore, an accelerated generalized minimax-concave (AGMC) sparse regularization method is proposed here.

ℓ₁ sparse regularization method

Because equation (6) has infinitely many solutions in the under-determined case, additional constraints are needed to get a satisfactory solution. Considering that the most entries in f are equal or close to zero, the most straightforward approach to obtain a sparse solution of equation (6) is to minimize the regularized linear least squares cost function with the ℓ₀-pseudo norm, mathematically expressed as

$\underset{\mathrm{f}}{\min }\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_0$

(7)

where ${\Vert \cdot \Vert }_0$ refers to the ℓ₀-pseudo norm counting the number of non-zero components. However, the optimal solution can be obtained by solving equation (7) via enumeration, which is an NP-hard problem. When the dimension of A becomes large, it is hard to obtain a satisfactory solution. Fortunately, replacing the ℓ₀-pseudo norm with the ℓ₁ norm in the regularization term can avoid this situation well while promoting the sparsity of the solution. Therefore, equation (7) can be relaxed as

$\underset{\mathrm{f}}{\min }\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_1$

(8)

where ${\Vert \mathbf{f}\Vert }_1={\displaystyle \sum }_{i=1}^n\left|f_i\right|$ denotes the ℓ₁ norm. Unlike equation (7), solving equation (8) is actually a convex optimization problem. At present, some efficient algorithms have been used to obtain the global optimal solution, including the interior point method, ¹¹ ISTA, ¹³ and so on. Algorithm 1 interprets the procedure of the well-known ISTA, where maxeig(A^TA) denotes the maximum eigenvalue of matrix A^TA, and soft(·) is defined as

$\text{soft}\left(\mathbf{f};q\right)=\text{sign}\left(\mathbf{f}\right)\cdot \max \left(\left|\mathbf{f}\right|-q,0\right)$

(9)

where f is the vector input and q is the threshold value. Although ℓ₁-norm minimization is feasible and widely used for inducing sparse solutions, this approach has been proved in theory and practice that it tends to underestimate the amplitude of the impact force.

Generalized minimax-concave sparse regularization method

In order to improve the precision and sparsity of the solution, an alternative non-convex regularizer has been developed in this study to replace the ℓ₁ norm. The cost function with the non-convex regularization term is generally expressed as

$\underset{\mathbf{f}}{\min }\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda \psi \left(\mathbf{f}\right)$

(10)

where ψ(f) is the non-convex penalty. Herein, we introduce a non-convex generalized minimax-concave (GMC) penalty into the impact force identification model in equation (10). First, a common non-convex penalty function ϕ _c named the scaled minimax-concave (MC) penalty ³³ is defined as

${\varphi }_c\left(f\right)=\{\begin{array}{ll}\left|f\right|-\frac{1}{2}{c}^2{f}^2,& \left|f\right|<1/c^2\\ \frac{1}{2c^2},& \left|f\right|\ge 1/c^2\end{array}$

(11)

Meanwhile, considering a univariate case, we define a function g as

$g\left(f\right)=\frac{1}{2}{\left(s-af\right)}^2+\lambda {\varphi }_c\left(f\right)$

(12)

where g(f) is convex if c² ≤ a²/λ (Proposition 3 in Ref. ²³ ). Then the minimum value of the convex function g(f) can be calculated by the firm threshold function ³⁴ as shown below:

$\text{firm}\left(f;\lambda ,\mu \right)=\{\begin{array}{ll}0,& \left|f\right|\le \lambda \\ \mu \left(\left|f\right|-\lambda \right)/\left(\mu -\lambda \right)\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(f\right),& \lambda \le \left|f\right|\le \mu \\ f,& \left|f\right|\ge \mu \end{array}$

(13)

where 0 < λ < μ. In order to better reflect the advantages of the MC penalty term, two diagrams are portrayed in Figure 1. One can find that the MC penalty is closer to the morphological features of the ℓ₀-psuedo norm penalty than the ℓ₁ norm. If f > λ, the value of the firm threshold function is closer to the exact value than the soft threshold function when searching for the minimum value of the cost function.

Figure 1.

One-dimensional examples of (a) penalties and (b) threshold functions.

Equivalently, equation (11) can be written in a Moreau envelope form as ³⁵

${\varphi }_c\left(f\right)=\left|f\right|-\underset{v\in \left\{0,f-1,f+1\right\}}{\min }\left\{\left|v\right|+\frac{1}{2c^2}{\left(f-v\right)}^2\right\}$

(14)

Next, extending the MC penalty in Equation (14) to a multivariate case, the GMC penalty is obtained, which is expressed as ²³

${\psi }_C\left(\mathbf{f}\right)={\Vert \mathbf{f}\Vert }_1-\underset{\mathbf{v}}{\min }\left\{{\Vert \mathbf{v}\Vert }_1+\frac{1}{2}{\Vert \mathbf{C}\left(\mathbf{f}-\mathbf{v}\right)\Vert }_2^2\right\}$

(15)

where the matrix C is employed to scale the penalty ψ _C . To visualize the GMC penalty, the matrix C is given as

$\mathbf{C}=\left[\begin{array}{ll}1& 0\\ 1/2& 1/2\\ 0& 1\end{array}\right]$

(16)

Then, the two-dimensional GMC penalty graph and its corresponding contour figure are shown in Figure 2. It can be found that the GMC penalty is a nonvex function.

Figure 2.

Two-dimensional examples of (a) the GMC penalty scaled by the matrix C and (b) the corresponding contour.

Substituting equation (15) into equation (10), one can get

$\begin{array}{l}G\left(\mathbf{f}\right)=\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\psi }_C\left(\mathbf{f}\right)\\ \begin{array}{l}=\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_1-\underset{\mathbf{v}}{\min }\left\{\lambda {\Vert \mathbf{v}\Vert }_1+\frac{\lambda }{2}{\Vert \mathbf{C}\left(\mathbf{f}-\mathbf{v}\right)\Vert }_2^2\right\}\\ =\underset{\mathbf{v}}{\max }\left\{\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_1-\lambda {\Vert \mathbf{v}\Vert }_1-\frac{\lambda }{2}{\Vert \mathbf{C}\left(\mathbf{f}-\mathbf{v}\right)\Vert }_2^2\right\}\\ \begin{array}{l}=\underset{\mathbf{v}}{\max }\left\{\frac{1}{2}{\mathbf{f}}^{\mathrm{T}}\left({\mathbf{A}}^{\mathrm{T}}\mathbf{A}-\lambda {\mathbf{C}}^{\mathrm{T}}\mathbf{C}\right)\mathbf{f}+\lambda {\Vert \mathbf{f}\Vert }_1+g\left(\mathbf{f},\mathbf{v}\right)\right\}\\ =\frac{1}{2}{\mathbf{f}}^{\mathrm{T}}\left({\mathbf{A}}^{\mathrm{T}}\mathbf{A}-\lambda {\mathbf{C}}^{\mathrm{T}}\mathbf{C}\right)\mathbf{f}+\lambda {\Vert \mathbf{f}\Vert }_1+\underset{\mathbf{v}}{\max }g\left(\mathbf{f},\mathbf{v}\right)\end{array}\end{array}\end{array}$

(17)

where g(f, v) is affine in f. The last term $\underset{\mathit{v}}{\max }g\left(\mathbf{f},\mathbf{v}\right)$ maintains convexity in f because it is the point-wise maximum of a set of convex functions indexed by v (Proposition 8.14 in Ref. ³⁵ ). Therefore, if the following condition equation (18) is satisfied, the cost function G(f) is convex (Theorem 1 in Ref. ²³ ). Then convex optimization algorithms can be sought to find the solution to equation (17).

${\mathbf{C}}^{\mathrm{T}}\mathbf{C}\underset{\_}{\prec }\frac{1}{\lambda }{\mathbf{A}}^{\mathrm{T}}\mathbf{A}$

(18)

In accordance with equation (18), once the matrix A has been obtained, it is possible to set the matrix C to a specific value as

$\mathbf{C}=\sqrt{\frac{\mathrm{\xi }}{\lambda }}\mathbf{A},0\le \mathrm{\xi }\le 1$

(19)

where ξ determines the non-sconvexity of the GMC penalty. If ξ = 0, the GMC penalty subsides to the ℓ₁ − norm penalty; if ξ = 1, the GMC penalty is maximally non-convex. Based on empirical observations, the value of ξ is typically selected from the interval $\left[\mathrm{0.5,0.8}\right]$ . ²³

Due to the absence of an analytical solution to equation (17), obtaining the global optimal solution requires the use of proximal algorithms, such as the forward-backward splitting (FBS) algorithm.^23,35 The convergence analysis has been presented in Refs.^23,36 Consequently, the problem of solving equation (17) can be reformulated as a saddle point problem which exhibits convexity with respect to f and concavity with respect to v, which can be expressed as

$\left({\mathbf{f}}^{\text{opt}},{\mathbf{v}}^{\text{opt}}\right)=\arg \underset{\mathbf{f}}{\min }\underset{\mathbf{v}}{\max }\left\{G\left(\mathbf{f},\mathbf{v}\right)\right\}$

(20)

where

$G\left(\mathbf{f},\mathbf{v}\right)=\frac{1}{2}{\Vert \mathbf{s}-\mathbf{A}\mathbf{f}\Vert }_2^2+\lambda {\Vert \mathbf{f}\Vert }_1-\lambda {\Vert \mathbf{v}\Vert }_1-\frac{\mathrm{\xi }}{2}{\Vert \mathbf{A}\left(\mathbf{f}-\mathbf{v}\right)\Vert }_2^2$

(21)

and the solution process of equation (20) via FBS is shown in Algorithm 2.

It is proved that the GMC penalty can effectively bypass the amplitude-underestimation problem compared with the ℓ₁-norm penalty. ²³ There exist two time-consuming gradient descent steps in Algorithm 2, especially when the dimension of the transfer matrix A is large. Therefore, an accelerated GMC sparse regularization method is proposed.

Accelerated generalized minimax-concave sparse regularization method

In light of the convex nature of the GMC regularized objective function, the Nesterov’s acceleration scheme is applied to the conventional GMC method to speed up convergence,^37,38 thereby improving the computational efficiency of the solution. The idea of this technique is to find a better point to calculate the update point. We are to cover some details of the acceleration process here. According to the soft threshold update rule, the updating of f^(k+1) and v^(k+1) is partially determined by w^(k) and u^(k), respectively. To expedite the convergence rate, we employ an extrapolation scheme to calculate w^(k+1) and u^(k+1) based on the present w^(k), u^(k) and their preceding values w^(k−1), u^(k−1). Accordingly, the calculation formulas can be expressed as

${\mathbf{w}}^{\left(k+1\right)}=\left(1-\delta \right){\mathbf{w}}^{\left(k\right)}+\delta {\mathbf{w}}^{\left(k-1\right)}$

(22)

${\mathbf{u}}^{\left(k+1\right)}=\left(1-\delta \right){\mathbf{u}}^{\left(k\right)}+\delta {\mathbf{u}}^{\left(k-1\right)}$

(23)

where

$\delta =\frac{1-d^{\left(k\right)}}{d^{\left(k+1\right)}},d^{\left(k+1\right)}=\frac{1+\sqrt{1+4{\left(d^{\left(k\right)}\right)}^2}}{2}$

Using this acceleration scheme, Algorithm 2 can be improved to Algorithm 3 as shown below. The main difference between Algorithm 3 and Algorithm 2 is that the iterative shrinkage operator soft(·) is employed at the extrapolation points w^(k+1) and u^(k+1) calculated by Equations (22) and (23) instead of the previous two points w^(k) and u^(k). As analyzed in Ref. ³⁷ , the complexity rate of Algorithm 3 is $O\left(1/k^2\right)$ , while that of Algorithm 2 is $O\left(1/k\right)$ .

K-sparsity criterion for the regularization parameter

It is often a difficult problem to choose an appropriate regularization parameter. Moreover, the choice of the parameter tends to directly affect the performance of the regularization method. There exist various criteria that can be used to choose the regularization parameter, such as Akaike Information Criterion (AIC), ³⁹ Bayesian Information Criterion (BIC), ⁴⁰ the L-curve method, ⁴¹ the generalized cross-validation method, ⁴² the K-sparsity method, ¹⁹ etc. Considering comprehensively, in order to take full advantage of sparse prior knowledge for impact forces and to also facilitate the uniform selection of regularization parameters, we adopt the K-sparsity criterion in the three algorithms mentioned above, namely, ISTA, GMC, and AGMC, for the adaptive selection of the regularization parameter λ.

It can be observed from the above three algorithms that λ only affects the threshold value in the soft thresholding function, and does not directly affect the iterative convergence speed. Given that λ determines how many coefficients are equal to zero in each iteration and how many should be preserved, we can make use of the K-sparsity prior knowledge of impact forces to adaptively set the value of λ at each iteration. The specific operation is to replace the parameter λμ with the adaptive threshold M^(k) which is set as the Kth largest coefficient in absolute value of the vector w^(k) in the ISTA and GMC algorithms, denoted as $M^{\left(k\right)}={\mathbf{w}}_{\left[K\right]}^{\left(k\right)}$ , and a similar substitution is also made in the AGMC algorithm, which is expressed as $M^{\left(k+1\right)}={\mathbf{w}}_{\left[K\right]}^{\left(k+1\right)}$ . Thus, Algorithm 1, Algorithm 2, and Algorithm 3 can be further updated to Algorithm 4, Algorithm 5, and Algorithm 6 as follows.

Numerical validation

In this section, to verify the performance of the proposed method AGMC in improving convergence speed and its reliability of impact force identification, a series of simulations are carried out. First, it is necessary to discuss the selection strategy of the regularization parameter K. Moreover, we compare the impact force identification results of AGMC, GMC, and ISTA when only changing the impact position with other conditions kept the same. Finally, we comparatively investigate the effect of noise on the accuracy of the three methods in reconstructing and localizing impact forces.

Problem description

To conduct numerical validation, a composite plate with clamped opposite edges is considered in this study. The plate has dimensions of 400mm × 300mm ×6mm and its material properties are provided in Table 1. Finite element model (FEM) in ANSYS is used to describe the structural dynamics of the plate subjected to external impact forces. The plate is discretized into 16 × 12 quadratic shell elements, and the first three natural frequencies are found to be 395.88 Hz, 441.90 Hz, and 621.25 Hz, respectively. To account for damping effects, the Rayleigh damping C = αM + βK is adopted with coefficients of α = 2.6240 and β = 3.7994 × 10⁻⁷.

Table 1.

Material properties of the composite plate.

Elastic modulus	Shear modulus	Poisson’s ratio	Density	Layer-ups
E1 = 135.0 GPa,	G12 = 4.47 GPa,	ν12 = 0.3,	ρ = 1560 kg/m3	$\left[45°/0°/-45°/90°\right]$
E2 = 8.8 GPa,	G23 = 3 GPa,	ν23 = 0.4,
E3 = 8.8 GPa	G13 = 4.47 GPa	ν13 = 0.3

The numerical validation is based on the strain measurements of the plate as structural responses. The arrangement of strain sensors and potential impact locations for all subsequent simulation studies is illustrated in Figure 3, where the optional measurement positions are labeled as S₁∼S₂₀, and the potential impact locations are numbered as P₁∼P₁₅. To simulate the under-determined cases, only four sensors are utilized to acquire the responses, as depicted in Figure 3.

Figure 3.

The composite plate with clamped opposite edges: (a) graphical representation of an applied impact force on the composite plate and (b) arrangement of strain gauges and potential impact locations. A total of 15 potential impact locations are taken into account, with only four selected measurement positions at S₈, S₁₀, S₁₄, and S₁₈.

A short-duration impulse force in the form of a Gaussian function is applied perpendicular to the plate, which can be expressed mathematically as

$f=B{e}^{-\frac{{\left(t-t_0\right)}^2}{2T^2}}$

(24)

where B denotes the amplitude of the impact force, t₀ represents the time at which the impact occurs, and T controls the duration of the impact. In the subsequent simulations, these parameters are set to t₀ = 0.0078 s, T = 9.7656e-5 s, and B = 150.

To generate synthetic signals that are representative of actual measured signals, Gaussian white noise with a standard deviation σ is added to the simulated dynamic responses, which is expressed as

$\stackrel{\sim }{\mathbf{s}}=\mathbf{s}+\sigma \mathbf{r}$

(25)

where the response vector $\stackrel{\sim }{\mathit{s}}$ represents the dynamic response corrupted by Gaussian white noise, s denotes the pure response vector, and the vector r is a sequence of independent random values drawn from a normal distribution with zero mean and unit standard deviation.

To obtain IRF h(t) in advance, we utilize the strain mode shapes of the plate and employ the Green function method. ⁴³ To assess the accuracy of the force reconstruction, the relative error (RE) between the actual force vector, denoted as f _p , and the estimated one ${\stackrel{\sim }{\mathbf{f}}}_p$ at the impact location p, can be defined as

$RE=\frac{{\Vert {\mathbf{f}}_p-{\stackrel{\sim }{\mathbf{f}}}_p\Vert }_2}{{\Vert {\mathbf{f}}_p\Vert }_2}\times 100\%$

(26)

where p denotes the sequential number of the specific location that is subjected to the impact force.

In addition, to serve as a local performance indicator for assessing the amplitude deviation in the identification results of impact forces, the peak relative error (PRE) is defined as

$PRE=\frac{{\Vert \max \left({\mathbf{f}}_p\right)-\max \left({\stackrel{\sim }{\mathbf{f}}}_p\right)\Vert }_2}{{\Vert \max \left({\mathbf{f}}_p\right)\Vert }_2}\times 100\%$

(27)

Furthermore, for the under-determined case, the localization error (LE) is required to assess the accuracy of the impact force localization, which is defined as

$LE=\left(1-\frac{{\Vert {\stackrel{\sim }{\mathbf{f}}}_p\Vert }_2}{{\displaystyle \sum }_{i=1}^{n_f}{\Vert {\stackrel{\sim }{\mathbf{f}}}_i\Vert }_2}\right)\times 100\%$

(28)

in which n _f is the total number of potential impact locations.

In order to facilitate comparison, it is necessary for the iteration termination criteria of the three algorithms to be consistent, which is denoted as

$\frac{{\Vert {\stackrel{\sim }{\mathbf{f}}}^{\left(k+1\right)}-{\stackrel{\sim }{\mathbf{f}}}^{\left(k\right)}\Vert }_1}{{\Vert {\stackrel{\sim }{\mathbf{f}}}^{\left(k\right)}\Vert }_1}\le \epsilon$

(29)

in which the value of ɛ is set to 1e-4.

Regularization parameter selection strategy

For choosing an appropriate regularization parameter K using the K-sparsity criterion, we discuss the effect of different K values on identification results of the three previously mentioned methods, namely, AGMC, GMC, and ISTA.

In this study, an impact force is applied perpendicularly to the plate at P₂. Moreover, parameters are set as the sampling frequency f _s = 10240Hz and the data analysis length N is 512 so that the size of transfer matrix A is 2048 × 7680 correspondingly. Gaussian white noise is added to all pure strain responses as an interference, and the signal-to-noise ratio (SNR) is set to 35, 30, 25, 20, 15, 10 dB, respectively. The K value varies from 1 to 50. In addition, the maximum iteration is fixed as 3000, and the non-convexity parameter ξ in both AGMC and GMC methods is fixed at 0.5. All programs are executed in MATLAB 2020a on a processor Inter^(R) Core^(R) i5-10400F, 24 GB RAM. To reduce the uncertainty influence, one hundred runs are conducted independently at each noise level.

Figures 4 and 5 show the relationship of RE values and LE values with respect to K values, respectively. It is discovered that REs of AGMC have been kept at the lowest and its localization accuracy remains the highest among the three methods when the K value varies from 6 to 50. Furthermore, AGMC can maintain its identification results at low REs and LEs over a large range of K values and this range is less affected by noise levels due to the great robustness of AGMC. It is worth noting that the RE variation trends of the three methods tend to be flat from K = 6 under 10 dB noise, and the rising trends of their LEs also slow down to a certain extent. The phenomenon of insensitivity to the value of K can also be observed in Ref. ¹⁹ According to the above phenomena, the K value of AGMC here can be chosen in the range of 6 to 20 under the requirement of low RE and LE. In the next simulation studies, the K value of the three methods is uniformly set as 15, which is almost consistent with the sparsity of the applied impact force in time domain.

Figure 4.

Relationship between RE and K values for identifying the same impact imposed at P₂ by AGMC, GMC, and ISTA at noise levels of 35, 30, 25, 20, 15, 10 dB, respectively.

Figure 5.

Relationship between LE and K values for identifying the same impact imposed at P₂ by AGMC, GMC, and ISTA at noise levels of 35, 30, 25, 20, 15, 10 dB, respectively.

Impact force identification using AGMC

To assess the performance of the proposed method in improving the convergence rate and identifying the impact force, we compare iterative convergence curves, force reconstruction, and localization accuracy of the proposed AGMC method with GMC and ISTA methods under the same simulation conditions. Utilizing the Monte Carlo method in this study, eight locations stochastically selected from the 15 potential impact locations are taken into consideration here. An impact force is applied to the eight locations, respectively; meanwhile, the noise level of measurement signals is set to 25 dB. Likewise, to remove randomness, one hundred runs are performed independently at each impact location.

Taken as examples, iterative convergence curves and identification results of the impact force applied at point P₃ and P₁₅, respectively, under the maximum iteration number of 3000 are shown in Figures 6 and 7. The cost value of AGMC increases significantly in the first iterations and then decreases rapidly. The final cost value of AGMC is very close to that of GMC. One can note that AGMC has the fastest convergence rate among the three methods. Furthermore, the amplitude and topography of AGMC identification results are closest to the exact forces, while the amplitude of the identification results obtained by ISTA differs the most from the true one. For further verification, relationships between the three accuracy indicators (PRE, RE, and LE) and maximum iteration number are portrayed in Figure 8. One can find that AGMC maintains the lowest PREs, REs, and LEs with the same maximum iteration number as GMC and ISTA and its three errors in the iterative process decline fastest among the three methods.

Figure 6.

Iterative convergence curves of AGMC, GMC, and ISTA with a maximum iteration number of 3000: (a) impact at P₃ and (b) impact at P₁₅.

Figure 7.

Reconstructed impact time histories at all monitored locations under 25 dB noise using four strain gages when (a) real impact acts on P₃ and (b) real impact acts on P₁₅.

Figure 8.

Relationship between the three accuracy indicators and maximum iteration number of AGMC, GMC, and ISTA: (a1)–(a3) impact at P₃ and (b1)–(b3) impact at P₁₅.

Figure 9 depicts the comparison of PREs, REs, and LEs in the identification results of the three methods at the eight randomly selected locations and Table 2 indicates the respective computing time and iteration number. It can be discovered that the iteration number of AGMC is without exception the least among the three methods and its computing time is also shortened compared to GMC. In addition, it is worth noting that ISTA exhibits PREs all above 10%, which can be attributed to the biased estimation caused by the ℓ₁-norm regularization method. Although the computing time of AGMC is longer than that of ISTA, AGMC can achieve impact force identification with high accuracy, and its PREs, REs, and LEs are less than those of GMC and ISTA.

Figure 9.

The three accuracy indicators of identification results of impact forces at eight different locations under 25 dB noise level by AGMC, GMC, and ISTA: (a) PRE, (b) RE, (c) LE.

Table 2.

Computing time and iterations of AGMC, GMC, and ISTA for identifying impact forces at different locations under 25 dB noise level.

Impact location i		P2	P6	P14	P12	P3	P8	P1	P15
Computing time (s)	AGMC	48.15	42.43	49.79	46.41	41.38	43.36	46.42	60.24
	GMC	67.30	55.36	69.38	58.45	45.50	50.01	60.36	76.47
	ISTA	21.81	19.93	20.45	17.51	16.06	19.93	20.81	23.84
Iteration number	AGMC	1872	1658	1860	1742	1618	1698	1838	1757
	GMC	2634	2172	2612	2223	1790	1978	2365	2228
	ISTA	2364	2172	2115	1882	1744	1980	2244	1941

The effect of noise level

To assess the robustness of the AGMC method against varying levels of noise, we analyze and compare the performance of the AGMC, GMC, and ISTA methods in constructing and localizing the impact force applied at P₁ and P₁₅ locations, respectively. Six different noise levels, namely, 35 dB, 30 dB, 25 dB, 20 dB, 15 dB, 10 dB, are considered in this study. Similarly, one hundred separate runs are conducted at each noise level to avoid randomness.

Table 3 represents the computing time and iteration number of identification results via AGMC, GMC, and ISTA at different noise levels. Figure 10 demonstrates the relationships between the three accuracy indicators (PRE, RE, and LE) and noise levels when impacting at the two locations, respectively. Undoubtedly, AGMC almost iterates minimally and it is more computationally efficient than GMC. Interestingly, AGMC still holds the highest reconstruction and localization accuracy at various noise levels. Additionally, tendencies in Figure 10 indicate that the reconstruction accuracy of AGMC decreases with the noise augmentation, while REs remain below 15% over a range of noise levels from 35 dB to 10 dB. The PEs of the AGMC identification results can generally be kept below 5% during the process of increasing noise, and only slightly increase when SNR is lower than 20 dB, but not more than 10%. Generally speaking, AGMC performs better than GMC and ISTA at various noise levels, in other words, it is quite robust to noise.

Table 3.

Computing time and iterations of identification results for impact forces imposed at P₁ and P₁₅ respectively by AGMC, GMC, and ISTA at different noise levels.

Position		P1						P15
SNR (dB)		35	30	25	20	15	10	35	30	25	20	15	10
Computing time (s)	AGMC	60.76	60.86	46.42	61.89	59.99	58.15	54.10	57.49	60.24	59.67	56.97	61.89
	GMC	80.07	80.53	60.36	80.11	82.65	82.04	70.81	74.94	76.47	74.64	77.89	77.93
	ISTA	26.81	26.94	20.81	26.99	26.39	24.22	23.09	23.47	23.84	23.32	23.17	22.66
Iteration number	AGMC	1785	1791	1838	1815	1759	1667	1725	1733	1757	1790	1725	1843
	GMC	2366	2363	2365	2369	2381	2346	2247	2243	2228	2254	2338	2330
	ISTA	2210	2214	2244	2175	2131	1975	1987	1969	1941	1923	1908	1899

Figure 10.

The three accuracy indicators of identification results for the same impact force via AGMC, GMC, and ISTA at different noise levels: (a1)–(a3) impact at P₁ and (b1)–(b3) impact at P₁₅.

Experimental verification

In this section, an experimental verification is conducted on a composite plate to validate the efficacy of the AGMC method in accelerating the convergence speed and improving the accuracy of impact force reconstruction and localization. An under-determined case is considered here, where only four strain gages are employed to monitor 15 potential impact locations, thus ensuring consistency in the dimension of the transfer matrix A with the simulation.

Experimental set-up

The experimental set-up involves an opposite-side-clamped composite plate, which possesses the same parameters and properties as the simulation model. The details of this configuration can be observed in Figure 11. Impact forces are generated using an impact hammer (PCB 086C03), and the impulsive signals are measured by a force sensor embedded in the hammer head. Strain responses are obtained via strain sensors (PCB 740B02). Both force and strain signals are acquired simultaneously using LMS SCADASIII Data Acquisition System with a sampling frequency of 10240 Hz. Potential impact positions are denoted by P₁∼P₁₅, and strain sensor locations are numbered as S₁∼S₁₀.

Figure 11.

Experimental configuration consisted of the composite plate with two opposite edges clamped. Four strain sensors, namely S₁, S₂, S₇, S₉, are selected to monitor the 15 potential impact locations.

In accordance with the simulation study, the strain responses measured at S₁, S₂, S₇, and S₉ are utilized to monitor all 15 potential impact positions. FRF a _ij (ω), representing the relationship between the output position i and the input location j, is obtained via impact testing using the LMS modal testing module. IRF a _ij (t), calculated by applying IFFT to a _ij (ω), is discretized to form the Toeplitz matrix. Consistent with the simulation, a data analysis length N = 512 is selected, leading to the dimension of A being 2048 × 7680. The ɛ value of the iteration termination criteria is set to 1e-4.

Results and discussion

Taking impact forces imposed at P₅ and P₁₁ as examples, the K value of the three methods is set as 15 and the non-convexity parameter ξ of AGMC and GMC is fixed as 0.5. Figures 12 and 13 show the time-history reconstruction and localization results of impact forces and iterative convergence curves of the three methods. Table 4 lists the PREs, REs, LEs, and computing time of the identification results. The iterative convergence curves in two figures intuitively illustrate that AGMC convergences faster than GMC and ISTA with fewer iterations. Table 4 also presents twice computing time of AGMC are 61.55 s and 76.54 s, which are 13.02% and 16.64% shorter than that of GMC when the impact occurs at P₅ and P₁₁, respectively.

Figure 12.

Identification results of the impact force applied to P₅ by AGMC, GMC and ISTA: (a)–(c) localization results of AGMC, GMC, and ISTA respectively, (d) time-history reconstruction results, and (e) iterative convergence curves of AGMC, GMC, and ISTA with 1802, 2024, and 2272 iterations, respectively.

Figure 13.

Identification results of the impact force applied to P₁₁ by AGMC, GMC and ISTA: (a)–(c) localization results of AGMC, GMC, and ISTA, respectively, (d) time-history reconstruction results, and (e) iterative convergence curves of AGMC, GMC, and ISTA with 2286, 2737, and 2498 iterations, respectively.

Table 4.

The three accuracy indicators and computing time of AGMC, GMC, and ISTA on identification results of impact forces applied to P₅ and P₁₁, respectively.

Position	PRE (%)			RE (%)			LE (%)			Computing time (s)
	AGMC	GMC	ISTA	AGMC	GMC	ISTA	AGMC	GMC	ISTA	AGMC	GMC	ISTA
P5	6.98	11.56	20.48	14.05	26.47	32.28	9.36	19.41	21.04	61.55	70.76	27.35
P11	4.24	9.77	20.98	11.01	23.54	32.41	8.05	18.67	21.45	76.54	91.82	29.32

Meanwhile, as shown in Figures 12 and 13, the force time-history reconstruction results of AGMC are the closest to the amplitude and shape of the exact forces and its force localization results are most accurate among the three methods. We noticed in Table 4 that AGMC holds the lowest PREs, REs, LEs, 6.98%, 14.05%, 9.36% at P₅, and 4.24%, 11.01%, 8.05% at P₁₁, respectively. It implies that when AGMC, GMC, and ISTA all reach the convergence criteria, AGMC can reconstruct impact forces with highest reconstruction and localization accuracy among the three methods with a shorter computing time than GMC. Meanwhile, although ISTA holds the shortest computational time, its REs, PREs, and LEs are significantly large. In particular, its PREs exceeding 20% are unacceptable identification results. This also highlights that ISTA is prone to the issue of underestimating the amplitude of the impact force in practical applications.

Figure 14 illustrates the relationships between the three accuracy indicators (PRE, RE, and LE) and maximum iteration number of AGMC, GMC, and ISTA when impacting at P₅ and P₁₁, respectively. One can note that under the same maximum iteration limit, AGMC consistently achieves the lowest values of RE, PRE, and LE among the three methods. Additionally, AGMC also exhibits the fastest convergence among the three methods with its PREs, REs, and LEs decreasing the steepest during the iteration process.

Figure 14.

Relationships between the three accuracy indicators and maximum iteration number of AGMC, GMC, and ISTA: (a1)–(a3) impact at P₅ and (b1)–(b3) impact at P₁₁.

From all the above results, we can conclude that the AGMC method effectively accelerates the convergence rate and thus solves the inefficiency problem of GMC. In addition, AGMC exhibits the highest accuracy among the three methods in terms of reconstructing and localizing impact forces under the under-determined condition, using the same iteration termination criteria.

Conclusions

In this paper, an accelerated generalized minimax-concave sparse regularization method is proposed for the purpose of reconstructing and localizing impact forces from limited structure responses under the under-determined condition. Although the GMC penalty is non-convex, it can still preserve the convexity of the cost function under a certain condition, so the global optimal solution can be obtained through convex optimization algorithms. In view of the fact that the solution efficiency of the GMC method decreases as the transfer matrix dimension increases, the Nesterov’s acceleration technique is utilized to improve the computational efficiency and the effectiveness of this strategy is verified systematically through simulations and experiments. The K-sparsity criterion is also used to adaptively select regularization parameters, and we noted that the K value can be selected in a relatively wide range. Simulations and experiments are carried out on a composite plate with two opposite edges clamped, both monitoring 15 potential impact locations with four strain sensors.

In summary, the convergence speed of AGMC is much faster than that of GMC in both simulations and experiments, which greatly improves the computational efficiency. Moreover, in the simulations, the AGMC method outperforms GMC and ISTA methods in force reconstruction and localization at different impact locations and different noise levels, with smallest PREs, REs, and LEs, even under 10 dB noise. Simulation results also demonstrate that AGMC is really resistant to noise. Finally, in the experimental study, AGMC still performs better than GMC and ISTA with highest reconstruction and localization accuracy. All these satisfactory performances of the AGMC method demonstrate its potential to solve the impact force identification problem under the under-determined condition in practical engineering applications.