CN105843780A - Sparse deconvolution method for impact load identification of mechanical structure - Google Patents

Abstract

The invention relates to a sparse deconvolution method for impact load identification of a mechanical structure. The method is used for solving the ill-posed nature of the impact load identification inverse problem and comprises steps as follows: 1) a frequency response function between an impact load acting point and a response measurement point of the mechanical structure is measured with a hammering method, a unit impulse response function is obtained through inverse fast fourier transform, discretization is further performed, and a transfer matrix is obtained; 2) an acceleration signal generated by impact load of the mechanical structure is measured with an acceleration sensor; 3) an L1-norm-based sparse deconvolution convex optimization model for impact load identification is established; 4) the sparse deconvolution optimization model is resolved with a primal-dual interior point method, and a sparse deconvolution solution, namely, a to-be-identified impact load, is obtained. The sparse deconvolution method is suitable for identifying the impact load acting on the mechanical structure. Compared with conventional Tikhonov regularization methods based on an L2 norm, the sparse deconvolution method has the advantages of high identification accuracy, high computation efficiency and high stability.

Description

Sparse deconvolution method for identifying impact load of mechanical structure

Technical Field

The invention belongs to the field of mechanical system state monitoring, and particularly relates to a sparse deconvolution method for identifying impact load of a mechanical structure.

Background

Impact loading is a very important type of dynamic loading, particularly in structural health monitoring of composite materials. For example, in the operation and maintenance process of a fan blade for wind power generation, the fan blade inevitably suffers from impact of wind sand, flying birds, hailstones, maintenance tools and other foreign objects, the occurrence frequency is high, and the impact damage is accumulated to cause potential safety hazards to the integrity and the bearing capacity of the fan blade; the impact of the sudden change airflow, the flying bird and the maintenance tool on the wings of the airplane. It is essential to monitor such impulse signals in real time, and it is very difficult to measure them directly.

The impact load identification method such as classical Tikhonov, truncated singular value decomposition and the like solves unknown load in an L2 norm space, relates to matrix inversion or decomposition operation, and is not suitable for large-scale inverse problem (the dimension of a transfer matrix is more than 10)⁴) The operation of (1); in recent years, function approximation methods (such as Db wavelets, Chebshiev polynomials, cubic B-splines, etc.) have been known according to the morphology of the impact loadThe multi-scholars are applied to single impact load identification of small-scale data volumes. The method achieves the aim of regularizing approximate dynamic load by controlling the number of basis functions, and the biggest dilemma of the application is that the regularization parameter, namely the number of the basis functions is difficult to determine, and the result is invalid due to too much or too little number of the basis functions; the neural network approach around impact load recognition requires a large number of samples to learn various types of impact events, which is very limited in practical engineering applications. The traditional regularization method is difficult to cope with continuous impact load identification under large data scale, and the condition number of a transfer matrix is very large and the numerical stability is very poor; in addition, when the response noise level is high, the performance of the traditional regularization method is reduced or even can not be improved.

Disclosure of Invention

Based on the above, the invention discloses a sparse deconvolution method for impact load identification of a mechanical structure, which comprises the following steps:

s100, measuring a frequency response function between an impact load action point of the mechanical structure and a response measuring point of the mechanical structure and calculating a transfer matrix;

s200, applying impact load to the mechanical structure and measuring impact response;

s300, constructing a sparse solution convolution convex optimization model of impact load identification based on the L1 norm based on the steps S100 and S200;

s400, solving the sparse deconvolution convex optimization model to obtain the sparse deconvolution solution of the impact load.

Compared with the prior art, the invention has the following advantages:

1. different from the traditional truncated singular value decomposition and Tikhonov regularization method based on the L2 norm, the sparse deconvolution method based on the L1 norm for impact load identification fully utilizes the time domain sparse characteristic of the impact load, and greatly inhibits the amplification of response noise in the identified impact load;

2. compared with the traditional Tikhonov regularization algorithm, the sparse deconvolution iterative algorithm has high identification precision, high calculation efficiency and strong stability;

3. the adopted primal-dual interior point method has the advantages that the inversion operation of the transfer matrix is not involved, and the regularization parameter is not required to be determined, and the iteration process is the regularization process;

4. the sparse deconvolution model and the corresponding primal-dual interior point method provided by the invention solve the problem of large data scale (the transfer matrix dimension is more than 10)⁴) And identifying the problem of the lower impact load.

Drawings

FIG. 1 is a flow chart of a sparse deconvolution method of impact load identification of a mechanical structure in one embodiment of the present invention;

FIG. 2 is a schematic view of a 300W composite wind turbine blade impact load identification apparatus according to one embodiment;

3(a), 3(b), 3(c) are 300W wind turbine blade data measured by force sensors and acceleration sensors in one embodiment, wherein FIG. 3(a) is the impact load signal, FIG. 3(b) is the acceleration response at point R1; FIG. 3(c) acceleration response at point R2;

4(a), 4(b) are graphs of the PDIPM identifying the trend of the dual gap variation with fan blade impact load for one embodiment, wherein the response signal of FIG. 4(a) is R1; FIG. 4(b) the response signal is R2;

FIGS. 5(a), 5(b) are sparse deconvolution of fan blade impact loads for one embodiment, where the response signal of FIG. 5(a) is R1; FIG. 5(b) the response signal is R2;

6(a), 6(b) are Tikhonov regularization results of fan blade impact loading in one embodiment, wherein FIG. 6(a) response signal is R1; the response signal of fig. 6(b) is R2.

Detailed Description

The invention will be further described with reference to fig. 1-6 and a specific embodiment, it being emphasized that the following description is merely exemplary and the invention is not limited in its application to the examples described below.

In one embodiment, the invention discloses a sparse deconvolution method of mechanical structure impact load identification, the method comprising the steps of:

s100, measuring a frequency response function between an impact load action point of the mechanical structure and a response measuring point of the mechanical structure and calculating a transfer matrix;

s200, applying impact load to the mechanical structure and measuring impact response;

s300, constructing a sparse solution convolution convex optimization model of impact load identification based on the L1 norm based on the steps S100 and S200;

s400, solving the sparse deconvolution convex optimization model to obtain the sparse deconvolution solution of the impact load.

In the present embodiment, the impact load may be applied to any position of the mechanical structure in the step S200, and is not limited to the measurement point in the step S100.

The method is used for solving the difficulties of low precision, low efficiency and poor stability of the traditional impact load identification method, and solves the impact load identification problem under the large data scale with high precision and high efficiency by adopting the advanced primal-dual interior point method in the field of convex optimization. According to the method, the sparse theory based on the L1 norm, which is widely concerned by the current scientific and engineering fields, is applied to the field of load identification, and the sparse convolution model is solved by adopting the advanced interior point method optimization algorithm in the optimization field, so that the amplification of response noise in the identified impact load is greatly inhibited.

In one embodiment, the step S100 specifically includes the following steps:

s1001, measuring a frequency response function between an impact load action point of a mechanical structure and a response measuring point of the mechanical structure;

s1002, performing fast Fourier inverse transformation on the frequency response function to obtain a unit impulse response function, and further performing discretization to obtain a transfer matrix.

In this embodiment, the method for measuring the frequency response function mainly includes a hammering method and a vibration exciter excitation method, where the hammering method is relatively convenient, and in this embodiment, the hammering method is preferably selected to measure the frequency response function.

In one embodiment, the step S200 is implemented by measuring an impact response generated by an impact load applied to the mechanical structure using a sensor.

In this embodiment, acceleration sensors are used to measure the corresponding signals generated by the dynamic loads acting on the mechanical structure, and velocity, displacement or strain sensors may also be used to measure the vibrational response.

In one embodiment, the sparse deconvolution convex optimization model in step S300 is:

$\begin{array}{cc}\underset{f}{\min imize}& \left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1\end{array};$

wherein | g | purple₂Representing the L2 norm of the vector, | | g | | non-volatile memory₁The L1 norm of the vector is represented, f represents the sparse deconvolution solution of the impact load, λ represents the regularization parameter, H system transfer matrix, and y is the impact load response vector.

In an embodiment, the sparse deconvolution convex optimization model in step S400 is solved by using a primal-dual interior point method (pdipmm), which specifically includes the following steps:

s401, constructing a dual gap function and a central path objective function by combining a sparse convolution convex optimization model, and initializing;

s402, calculating the searching direction of the central path target function;

s403, determining an iteration step length of the central path function;

s404, updating obstacle parameters in the central path function;

s405, updating the current solution based on the steps S401-404;

s406, setting a PDIPM iteration termination criterion, and if the current solution meets the iteration termination criterion, terminating the iteration process, wherein the current solution is a sparse deconvolution solution f of the impact load; otherwise, the iterative process returns to step S401 to continue the iterative computation until the pdipmm iteration termination criterion is satisfied.

In this embodiment, the dual gap function is used as a condition for determining the termination of the iteration of the dual interior point method, and the central path objective function is an objective function to be optimized by the dual interior point method.

In one embodiment, the dual gap function in step S401 is:

wherein the dual feasible variable v ═ 2(Hf-y), wherein | | | g | | survival₂Representing the L2 norm of the vector, | | g | | non-volatile memory₁L1 norm representing vector, f is sparse convolution solution, λ represents regularization parameter, H is system transfer matrix, y is load response vector, G (v) represents Lagrange dual problem, $\begin{array}{cc}\underset{v}{\max imize}G\left(v\right)=-\frac{1}{4}{v}^{T}v-v^{T}y;subjectto& {\left(A^{T}v\right)}_i\≤\mathrm{\λ}\end{array};$

the center path function is:

$\begin{array}{cc}\underset{(f,u)}{\min imize}& {\mathrm{\φ}}_t(f,u)=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}\underset{i=1}{\overset{n}{\Σ}}{u}_i+\frac{1}{t}\mathrm{\Φ}(f,u)\end{array}$

wherein t ∈ (0, ∞) is barrier parameter, u ∈ RⁿFor barrier constrained variable, u_iConstraining variable u ∈ R for a barrierⁿPhi (f, u) is a logarithmic barrier function;

said step (c) isThe initialization in S401 specifically includes: initializing a non-negative initial solution f ═ 0]^T∈RⁿBarrier constraint variable u ═ 1]^T∈RⁿThe termination threshold is more than 0, and the regularization parameter lambda is 0.02H^Ty||_∞And the initial obstacle parameter t is 1/lambda.

In one embodiment, the step S402 calculates the search direction [ Δ f ] of the center path function by using the following formula^T，Δu^T]^T：

${\▿}^2{\mathrm{\φ}}_t(f,u)\left[\begin{array}{c}\mathrm{\Δ}f\\ \mathrm{\Δ}u\end{array}\right]=\▿{\mathrm{\φ}}_t(f,u)$

Wherein,is a center path function phi_tA Hessian matrix of (f, u),is a center path function phi_tA gradient of (f, u);

setting the iteration step size in step S403 to S, and making S equal to β^jWherein j represents the smallest integer satisfying the following formula:

${\mathrm{\φ}}_t(f+{\mathrm{\β}}^j,u+{\mathrm{\β}}^j)\≤{\mathrm{\φ}}_t(f,u)+{\mathrm{\α\β}}^j\▿{\mathrm{\φ}}_t{(f,u)}^T{\[{\mathrm{\Δf}}^T,{\mathrm{\Δu}}^T\]}^T$

where α ∈ (0, 1/2) and β ∈ (0, 1) are straight line search parameters.

In the present embodiment, more preferable parameters are set as: α ═ 0.01 and β ═ 0.5.

In one embodiment, the step S404 is to update the barrier parameter t by using the following formula:

$t=\left\{\begin{array}{cc}max\left\{min\right\{2\mathrm{\μ}n/\mathrm{\η},\mathrm{\μ}t\},t\},& s\≥s_{min}\\ t,& otherwise\end{array}\right.$

where n is the dimension of the transfer matrix H and the constant μ 2, s_min＝0.5。

In one embodiment, the current solution in S405Including a current impact load sparse deconvolution solution f_newAnd the current solution u of the barrier constraint variable_newThe specific updating formula is as follows:

f_new＝f_old+sΔf

u_new＝u_old+sΔu

wherein f is_oldAnd u_oldThe last PDIPM iteration result.

In one embodiment, the pdipmm iteration termination criteria in step S406 are:

wherein the termination threshold represents an acceptable error.

In this embodiment, the value of the termination threshold is continuously tried, and if the value is too small, although the precision is slightly improved, the iteration step number is more and the calculation time is longer; conversely, a too large value may result in premature termination of the iteration; typical values are e.g. ═ 0.01.

In one embodiment, the invention discloses a sparse deconvolution method for impact load identification of a mechanical structure, which specifically comprises the following steps:

1) measuring the frequency response function and calculating the transfer matrix. Measuring a frequency response function H (omega) between an impact load action point of a mechanical structure and a response measuring point of the mechanical structure by adopting a hammering method, obtaining a unit impulse response function H (t) through Inverse Fast Fourier Transform (IFFT), and further discretizing to obtain a transfer matrix H, wherein omega represents a circular frequency variable, and t represents a time variable;

2) applying an impact load and measuring an acceleration impact response, and measuring a response signal y generated by the impact load acting on the mechanical structure using an acceleration sensor. Note that: the acceleration sensors in the step 1) and the step 2) can be replaced by speed, displacement and strain gauge sensors;

3) considering that the impact load signal is a sparse sequence in the time domain (in the test time, the impact force has a larger value only in the impact loading stage, and the non-loading region can be regarded as zero), constructing a sparse solution convolution convex optimization model based on the L1 norm:

$\begin{array}{cc}\underset{f}{\min imize}& \left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1\end{array};---\left(1\right)$

wherein | g | purple₂Representing the L2 norm of the vector, | | g | | non-volatile memory₁An L2 norm representing the vector, f representing the impact load, and λ representing the regularization parameter; unlike the regularization model energy minimization criterion based on the L2 norm, the regularization model based on the L1 norm suppresses the diffusion of the learned energy at all time points, preserving the sparsity characteristics of the signal to be reconstructed.

4) Compared with the traditional Tikhonov method based on L2 regularization, the method has the advantages of explicit solution f ═ (H)^TH+λI)^-1H^Ty is different, the L1 regularization model has no explicit solution, and a primal-dual interior point method is used for solving the sparse deconvolution optimization model, which specifically comprises the following steps:

initialization 4.1): combining the sparse deconvolution model in the formula (1), constructing a dual gap function:

$\mathrm{\η}=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1-G\left(v\right)---\left(2\right)$

wherein the dual feasible variable v ═ 2(Af-y) and Lagrange dual problem g (v):

$\begin{array}{cc}\underset{v}{\max imize}G\left(v\right)=-\frac{1}{4}{v}^{T}v-v^{T}y;subjectto& {\left(A^{T}v\right)}_i\≤\mathrm{\λ}\end{array}---\left(3\right)$

initialization 4.2): constructing a center path function:

$\begin{array}{cc}\underset{(f,u)}{\min imize}& {\mathrm{\φ}}_t(f,u)=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}\underset{i=1}{\overset{n}{\Σ}}{u}_i+\frac{1}{t}\mathrm{\Φ}(f,u)\end{array}---\left(4\right)$

wherein t ∈ (0, ∞) is barrier parameter, u ∈ RⁿConstraint variables for the obstacle and the logarithmic obstacle function Φ (f, u):

$\mathrm{\Φ}(f,u)=-\underset{i=1}{\overset{n}{\Σ}}ln(u_i+f_i)-\underset{i=1}{\overset{n}{\Σ}}ln(u_i-f_i)---\left(5\right)$

initialization 4.3): non-negative initial solution f ═ 0,.. 0, 0]^T∈RⁿBarrier constraint variable u ═ 1]^T∈RⁿThe termination threshold is more than 0, and the regularization parameter lambda is 0.02H^Ty||_∞The initial obstacle parameter t is 1/lambda; wherein | · | purple sweet_∞Represents an infinite norm;

step 4.1): calculating the search direction [ Delta f ] of the center path function^T，Δu^T]^T：

${\▿}^2{\mathrm{\φ}}_t(f,u)\left[\begin{array}{c}\mathrm{\Δ}f\\ \mathrm{\Δ}u\end{array}\right]=\▿{\mathrm{\φ}}_t(f,u)---\left(6\right)$

Wherein,is the Hessian matrix of the center path function,is the gradient of the central path function.

Step 4.2) calculating an iteration step s, and enabling s to be β^jWherein j is the smallest integer satisfying the following formula:

${\mathrm{\φ}}_t(f+{\mathrm{\β}}^j,u+{\mathrm{\β}}^j)\≤{\mathrm{\φ}}_t(f,u)+{\mathrm{\α\β}}^j\▿{\mathrm{\φ}}_t{(f,u)}^T{\[{\mathrm{\Δf}}^T,{\mathrm{\Δu}}^T\]}^T---\left(7\right)$

where α ∈ (0, 1/2) and β ∈ (0, 1) are straight line search parameters, typically with values α ═ 0.01 and β ═ 0.5.

Step 4.3): updating obstacle parameters:

$t=\left\{\begin{array}{cc}max\left\{min\right\{2\mathrm{\μ}n/\mathrm{\η},\mathrm{\μ}t\},t\},& s\≥s_{min}\\ t,& otherwise\end{array}\right.---\left(8\right)$

where n is the dimension of the transfer matrix H, the constants μ 2 and s_min＝0.5。

Step 4.4): and updating the current solution:

$\begin{array}{c}{f}_{new}=f_{old}+s\mathrm{\Δ}f\\ u_{new}=u_{old}+s\mathrm{\Δ}u\end{array}---\left(9\right)$

step 4.5): the ratio of the dual gap function and the dual objective function is taken as the PDIPM iteration termination criterion:

$\mathrm{\ϵ}\≤\frac{\mathrm{\η}}{G\left(v\right)}---\left(10\right)$

where the termination threshold represents an acceptable error, typically 0.01.

If the current solution f_newIf the above formula iteration termination criterion is met, terminating the iteration process to obtain a sparse deconvolution solution f of the impact load; otherwise, the iterative process returns to step 4.1) to continue the iterative computation until the above formula is satisfied.

FIG. 1 is a flow chart of a sparse deconvolution method for mechanical structure impact load identification, which is implemented by the present invention, and the method makes full use of the time domain sparsity of the impact load, constructs an impact load sparse deconvolution model according to the measured frequency response function and acceleration signal, and solves the problem by using an advanced primal-dual interior point method in the optimization field to realize the impact load identification, and comprises the following specific steps:

1) measuring the frequency response function and calculating the transfer matrix. Measuring a frequency response function H (omega) between a dynamic load application position point of a mechanical structure and a mechanical structure response measuring point by adopting a hammering method, obtaining a unit impulse response function H (t) through Inverse Fast Fourier Transform (IFFT), and further discretizing to obtain a transfer matrix H, wherein omega represents a circular frequency variable, and t represents a time variable; wherein, the hammering method is a test mode test method commonly used in the field;

11) test forThe experimental diagram is shown in FIG. 2, and the selected composite material fan blade has a length of 702mm, a maximum width of 112mm and a thickness of 11 mm. In order to simulate the real installation condition of the fan blade, the fixed support section of the cantilever beam support is used for clamping the position of 0-60 mm of the root of the fan blade, and two acceleration sensors with the models of PCB 333B50 are respectively arranged at the positions 120-160 mm away from the root and the end and are respectively marked as measuring points R1 and R2. The location of the simulated impact force is 70% from the blade root. The impact point is repeatedly knocked five times by adopting a pulse force hammer (a force sensor is embedded at the top of a hammer head) with the model of PCB 086C02, an impact force and an acceleration signal are synchronously recorded by a LMSSCADASIII data acquisition system, and a frequency response function from the impact load action point five times to an acceleration measuring point is H₁(ω)、H₂(ω)、H₃(ω)、H₄(omega) and H₅(ω), the average value of which is H (ω) calculated by the module LMS IMPACT;

12) the sampling frequency when measuring the system frequency response function is 2048Hz, the sampling time is 1s, and the data length is 2050. The condition numbers of the high dimensional transfer matrices for excitation point F and response points R1 and R2 are 4.92E +17 and 2.66E +18, respectively (the condition numbers are an index for measuring the degree of matrix morbidity). It is known that the fan blade impact load identification problem is severely ill.

2) Applying an impact load and measuring acceleration impact response, and measuring a response signal y generated by the impact load acting on the mechanical structure by using an acceleration sensor;

21) the fan blade susceptible impact position is tapped five times continuously by using an impact force hammer, and an acceleration signal and an impact load signal are synchronously recorded at a sampling frequency of 2048Hz by LMS SCADASIII data acquisition (as shown in figures 3(a) and 3(b) and 3 (c)). And the actually measured force signal is used as a comparison object of the identification result of the impact load sparse convolution method. Note that the point of action of this step to apply the impact is consistent with the point of action of the step 1) to measure the frequency response function, while the acceleration position remains constant throughout.

22) In the impact test, the impact hammer strikes the fan blade five times in succession, and is marked as F1, F2, F3, F4 and F5 respectively, as shown in FIG. 3 (a). To verify the efficiency of the sparse deconvolution method of the present invention, four consecutive impact loads in fig. 3(a) would be identified simultaneously. The duration of the continuous impact was 14s, the sampling frequency was 2050Hz, and the response data length was 26638.

23) Because the continuous impact time history is reconstructed at the same time, the dimension of the solved problem becomes very large, the dimension of the original transfer function is far smaller than the length of the actually measured response data, and the load identification control equation at the moment is defined as:

$\left[\begin{array}{c}y\left(\mathrm{\Δ}t\right)\\ y\left(2\mathrm{\Δ}t\right)\\ M\\ y\left(n\mathrm{\Δ}t\right)\\ y\left(\right(n+1\left)\mathrm{\Δ}t\right)\\ M\\ y\left(m\mathrm{\Δ}t\right)\end{array}\right]=\mathrm{\Δ}t\left[\begin{array}{ccccc}h\left(\mathrm{\Δ}t\right)& 0& L& 0& 0\\ h\left(2\mathrm{\Δ}t\right)& h\left(\mathrm{\Δ}t\right)& L& 0& 0\\ M& M& L& M& M\\ h\left(n\mathrm{\Δ}t\right)& h\left(\right(n-1)\mathrm{\Δ}t)& L& 0& 0\\ 0& h\left(n\mathrm{\Δ}t\right)& L& 0& 0\\ M& M& L& h\left(\mathrm{\Δ}t\right)& 0\\ 0& 0& L& h\left(2\mathrm{\Δ}t\right)& h\left(\mathrm{\Δ}t\right)\end{array}\right]\left[\begin{array}{c}f\left(\mathrm{\Δ}t\right)\\ f\left(2\mathrm{\Δ}t\right)\\ M\\ f\left(n\mathrm{\Δ}t\right)\\ f\left(\right(n+1\left)\mathrm{\Δ}t\right)\\ M\\ f\left(m\mathrm{\Δ}t\right)\end{array}\right]---\left(1\right)$

wherein: m is the data length of the measured response, n is the unit impulse response function, and m > n. After discretization, the deconvolution model of equation (1) can be rewritten into a compact form of a matrix-vector:

y＝Hf； (2)

wherein the new transfer matrix H ∈ R^m×mIn the sparse framework, the transfer matrix is called the observation matrix H. in this case, the response data dimension m for inverting dynamic loads is 26638, and the observation matrix H dimension is 26638 × 26638.

3) Considering that the impact load signal is a sparse sequence in the time domain (in the test time, the impact force has a larger value only in the impact loading stage, and the non-loading region can be regarded as zero), constructing a sparse solution convolution convex optimization model based on the L1 norm:

$\begin{array}{cc}\underset{f}{\min imize}& \left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1\end{array};---\left(4\right)$

wherein | g | purple₂Representing the L2 norm of the vector, | | g | | non-volatile memory₁An L1 norm representing the vector, f represents a sparse deconvolution solution of the impact load, and λ represents a regularization parameter;

4) compared with the traditional Tikhonov method based on L2 regularization, the method has the advantages of explicit solution f ═ (H)^TH+λI)^-1H^Ty is different, the L1 regularization model has no explicit solution, and a primal-dual interior point method is used for solving the sparse deconvolution optimization model, which specifically comprises the following steps:

initialization 4.1): in combination with the sparse deconvolution model of claim 1, a dual gap function is constructed:

$\mathrm{\η}=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1-G\left(v\right)---\left(5\right)$

wherein the dual feasible variable v ═ 2(Af-y) and Lagrange dual problem g (v):

$\begin{array}{cc}\underset{v}{\max imize}G\left(v\right)=-\frac{1}{4}{v}^{T}v-v^{T}y;subjectto& {\left(A^{T}v\right)}_i\≤\mathrm{\λ}\end{array}---\left(6\right)$

initialization 4.2): constructing a center path function:

$\begin{array}{cc}\underset{(f,u)}{\min imize}& {\mathrm{\φ}}_t(f,u)=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}\underset{i=1}{\overset{n}{\Σ}}{u}_i+\frac{1}{t}\mathrm{\Φ}(f,u)\end{array}---\left(7\right)$

wherein t ∈ (0, ∞) is an obstacle parameter, u ∈ RⁿConstraint variables for the obstacle and the logarithmic obstacle function Φ (f, u):

$\mathrm{\Φ}(f,u)=-\underset{i=1}{\overset{n}{\Σ}}\ln (u_i+f_i)-\underset{i=1}{\overset{n}{\Σ}}\ln (u_i-f_i)---\left(8\right)$

initialization 4.3): non-negative initial solution f ═ 0,.. 0, 0]^T∈RⁿBarrier constraint variable u ═ 1]^T∈RⁿThe termination threshold is more than 0, and the regularization parameter lambda is 0.02H^Ty||_∞The initial obstacle parameter t is 1/lambda; therein,. mu.g|·||_∞Represents an infinite norm; .

In the case, the feasibility of calculating the high-dimensionality inverse problem by the Tikhonov method is considered, the computing environment is a workstation which is configured into 2 Intel Xeon E5580 CPUs and 1 memory of 48 GB;

step 4.1): calculating the search direction [ Delta f ] of the center path function^T，Δu^T]^T：

${\▿}^2{\mathrm{\φ}}_t(f,u)\left[\begin{array}{c}\mathrm{\Δ}f\\ \mathrm{\Δ}u\end{array}\right]=\▿{\mathrm{\φ}}_t(f,u)---\left(9\right)$

Wherein,is the Hessian matrix of the center path function,is the gradient of the central path function.

Step 4.2) calculating an iteration step s, and enabling s to be β^jWherein j is the smallest integer satisfying the following formula:

${\mathrm{\φ}}_t(f+{\mathrm{\β}}^j,u+{\mathrm{\β}}^j)\≤{\mathrm{\φ}}_t(f,u)+{\mathrm{\α\β}}^j\▿{\mathrm{\φ}}_t{(f,u)}^T{\[{\mathrm{\Δf}}^T,{\mathrm{\Δu}}^T\]}^T---\left(10\right)$

where α ∈ (0, 1/2) and β ∈ (0, 1) are straight line search parameters, typically with values α ═ 0.01 and β ═ 0.5.

Step 4.3): updating obstacle parameters:

$t=\left\{\begin{array}{cc}max\left\{min\right\{2\mathrm{\μ}n/\mathrm{\η},\mathrm{\μ}t\},t\},& s\≥s_{min}\\ t,& otherwise\end{array}\right.---\left(11\right)$

where n is the dimension of the transfer matrix H, the constants μ 2 and s_min＝0.5。

Step 4.4): the updated current solution includes the current of the impact loadSolution f_newAnd the current solution u of the barrier constraint variable_new：

$\begin{array}{c}{f}_{new}=f_{old}+s\mathrm{\Δ}f\\ u_{new}=u_{old}+s\mathrm{\Δ}u\end{array}---\left(12\right)$

The current solution includes a current solution f of the impact load_newAnd the current solution u of the barrier constraint variable_new。

Step 4.5): the ratio of the dual gap function and the dual objective function is taken as the PDIPM iteration termination criterion:

$\mathrm{\ϵ}\≤\frac{\mathrm{\η}}{G\left(v\right)}---\left(13\right)$

where the termination threshold represents an acceptable error, typically 0.01.

If the current solution f_newIf the above formula iteration termination criterion is met, terminating the iteration process to obtain a sparse deconvolution solution f of the impact load; otherwise, the iterative process returns to step 4.1) to continue the iterative computation until the above formula is satisfied.

5) In order to quantitatively evaluate the accuracy of loads identified by sparse deconvolution algorithms PDIPM and a Tikhonov regularization method, respectively defining time domain global relative errors and impact load peak value relative errors as follows:

$\frac{\left|\right|f_{measured}-f_{identified}|{|}_2}{\left|\right|f_{measured}|{|}_2}\×100\%---\left(14\right)$

$\frac{\left|\right|\max \left(f_{measured}\right)-\max \left(f_{identified}\right)|{|}_2}{\left|\right|\max \left(f_{measured}\right)|{|}_2}\×100\%---\left(15\right)$

wherein f is_measuredAnd f_identifiedThe measured impact load of the force sensor and the impact load reconstructed by applying the regularization method are respectively.

51) Fig. 4(a), 4(b) show that when the sparse deconvolution algorithm PDIPM identifies successive impact loads acting on the thin shell structure using response signals R1 and R2, respectively, the dual gap decreases rapidly and then gradually flattens as the number of iteration steps increases. The number of terminating iteration steps is 23 and 26, respectively.

52) Fig. 5(a), 5(b) are the continuous impact loads reconstructed by the sparse deconvolution algorithm pdippm using the response signals R1 and R2. It can be seen that the five impact loads reconstructed by the PDIPM using the response signals R1 and R2, respectively, are highly matched with the measured load, and the reconstructed impact loads are very sparse in the impact duration interval. Table 1 gives that the relative error of the continuous impact load reconstructed by the PDIPM with the response signals R1 and R2 is only 13.02% and 12.93%, respectively. The peak forces of the five measured impact loads were 11.19N, 36.24N, 23.88N, 38.90N, and 26.00N, respectively. The peak forces of the load inverted with the response signal R1 (see fig. 5(a)) were 11.43N, 36.81N, 24.61N, 39.22N, and 26.20N, respectively, and the relative errors were 2.14%, 1.51%, 3.06%, 0.82%, and 0.77%, respectively (see table 1).

53) Fig. 6(a), fig. (b) are the continuous impact loads reconstructed by the Tikhonov regularization method with response signals R1 and R2. It can be seen that the Tikhonov solution is amplified in the unloaded region of the impact load, particularly when the response point R2 is used to invert the dynamic load (see fig. 6 (b)). Table 1 shows that the relative error of the continuous impact load reconstructed by the Tikhonov method with response signals R1 and R2 is as high as 58.09% and 290.24%, respectively. In contrast, with the aid of the prior information of the temporal sparsity of the impact load, the continuous impact load reconstructed by the PDIPM is almost zero in the impact unloaded region (see fig. 5(a), 5 (b)). It is worth noting that: in this case, although the relative error of the Tikhonov method is large, the peak relative error is still acceptable (less than 10%) or even smaller than the PDIPM peak relative error. As for the time required for the continuous impact load reconstruction, Tikhonov takes much more time than pdippm. Therefore, the sparse deconvolution method can stably solve the inverse problem of high-dimensionality and serious-morbidity load identification with high precision and high efficiency.

Table 1 comparison of PDIPM and Tikhonov identification results for continuous impact loads of fan blades

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents and improvements made within the spirit and principle of the present invention are intended to be included within the scope of the present invention.

Claims (10)

1. A sparse deconvolution method for impact load identification of a mechanical structure, the method comprising the steps of:

s100, measuring a frequency response function between an impact load action point of the mechanical structure and a response measuring point of the mechanical structure and calculating a transfer matrix;

s200, applying impact load to the mechanical structure and measuring impact response;

s300, constructing a sparse solution convolution convex optimization model of impact load identification based on the L1 norm based on the steps S100 and S200;

s400, solving the sparse deconvolution convex optimization model to obtain the sparse deconvolution solution of the impact load to be identified.

2. The method according to claim 1, wherein step S100 preferably comprises the following steps:

s1001, measuring a frequency response function between an impact load action point of a mechanical structure and a response measuring point of the mechanical structure;

s1002, performing fast Fourier inverse transformation on the frequency response function to obtain a unit impulse response function, and further performing discretization to obtain a transfer matrix.

3. The method of claim 1, wherein: in step S200, an impact response generated by an impact load applied to the mechanical structure is measured by a sensor.

4. The method of claim 1, wherein the sparse deconvolution convex optimization model in step S300 is:

$\begin{array}{cc}\underset{f}{\min imize}& \left|\right|Hf-y|{|}_2^2+\mathrm{\λ}|\left|f\right|{|}_1\end{array};$

wherein | g | purple₂Representing the L2 norm of the vector, | | g | | non-volatile memory₁And L1 norm representing the vector, f represents the impact load sparse deconvolution, lambda represents the regularization parameter, H is the transfer matrix, and y is the impact load response vector.

5. The method according to claim 4, wherein the step S400 solves the sparse deconvolution convex optimization model by using a primal-dual interior point method (PDIPM) to obtain a sparse deconvolution solution of the impact load, and specifically comprises the steps of:

s401, constructing a dual gap function and a central path objective function by combining a sparse convolution convex optimization model, and initializing;

s402, calculating the searching direction of the central path target function;

s403, calculating the iteration step length of the central path function;

s404, updating obstacle parameters in the central path function;

s405, iteratively updating the current solution based on the steps S401-404;

s406, setting a PDIPM iteration termination criterion, and if the current solution meets the iteration termination criterion, terminating the iteration process, wherein the current solution is a sparse deconvolution solution of the impact load f; otherwise, the iterative process returns to step S401 to continue the iterative computation until the pdipmm iteration termination criterion is satisfied.

6. The method of claim 5,

the dual gap function in step S401 is:

wherein, the dual feasible variable v is 2(Hf-y), f is the impact load sparse deconvolution, λ represents the regularization parameter, H is the system transfer matrix, y is the impact load response vector, G (v) represents Lagrange dual function,

$\begin{array}{ccc}\underset{v}{\max imize}& G\left(v\right)=-\frac{1}{4}{v}^{T}v-v^{T}y;subjectto& {\left(A^{T}v\right)}_i\≤\mathrm{\λ}\end{array};$

the center path objective function is:

$\begin{array}{cc}\underset{(f,u)}{\min imize}& {\mathrm{\φ}}_t(f,u)=\left|\right|Hf-y|{|}_2^2+\mathrm{\λ}\underset{i=1}{\overset{n}{\Σ}}{u}_i+\frac{1}{t}\mathrm{\Φ}(f,u)\end{array}$

wherein t ∈ (0, ∞) is barrier parameter, u ∈ RⁿFor barrier constrained variable, u_iConstraining variable u ∈ R for a barrierⁿPhi (f, u) is a logarithmic barrier function;

the initialization in step S401 specifically includes: initializing impact load sparse deconvolution solution f ═ 0]^T∈RⁿBarrier constraint variable u ═ 1]^T∈RⁿThe termination threshold is more than 0, and the regularization parameter lambda is 0.02H^Ty||_∞And the initial obstacle parameter t is 1/lambda.

7. The method of claim 6,

in step S402, the following formula is used to calculate the search direction [ Δ f ] of the center path function^T，Δu^T]^T：

${\▿}^2{\mathrm{\φ}}_t(f,u)\left[\begin{array}{c}\mathrm{\Δ}f\\ \mathrm{\Δ}u\end{array}\right]=\▿{\mathrm{\φ}}_t(f,u)$

Wherein,is a center path function phi_tA Hessian matrix of (f, u),is a center path function phi_tA gradient of (f, u);

the iteration step size in step S403 is S, where S is β^jWherein j represents the smallest integer satisfying the following formula:

${\mathrm{\φ}}_t(f+{\mathrm{\β}}^j,u+{\mathrm{\β}}^j)\≤{\mathrm{\φ}}_t(f,u)+{\mathrm{\α\β}}^j\▿{\mathrm{\φ}}_t{(f,u)}^T{\[{\mathrm{\Δf}}^T,{\mathrm{\Δu}}^T\]}^T$

where α ∈ (0, 1/2) and β ∈ (0, 1) are straight line search parameters.

8. The method of claim 7,

in step S404, the obstacle parameter t is updated by using the following formula:

$t=\left\{\begin{array}{cc}max\left\{min\right\{2\mathrm{\μ}n/\mathrm{\η},\mathrm{\μ}t\},t\},& s\≥s_{min}\\ t,& otherwise\end{array}\right.$

where n is the dimension of the transfer matrix H, with the constant μ 2, S_min＝0.5。

9. The method of claim 8,

the current solution in S405 comprises a current impact load sparse convolution solution f_newAnd the current solution u of the barrier constraint variable_newThe specific updating formula is as follows:

f_new＝f_old+sΔf

u_new＝u_old+sΔu

wherein f is_oldAnd u_oldDeconvoluting convolution solution and obstacle constraint for impact load sparsity in last PDIPM iteration resultThe result of the variable.

10. The method of claim 9,

the PDIPM iteration termination criterion in step S406 is: