A method for simultaneously identifying structural dynamic load and local damage by parallel computation

Abstract

This invention discloses a parallel computational method for simultaneous identification of structural dynamic loads and local damage, comprising the following steps: 1) Arranging accelerometers on the surface of a beam structure to collect structural responses, and storing the structural responses in a two-dimensional matrix form using a moving time window; 2) Establishing a finite element model using beam elements, and dividing the beam elements into neg groups to construct neg independently runnable identification subtasks; 3) For each subtask, using the element damage factor contained in that group as an optimization parameter, and combining matrix regularization to construct an optimization objective function; 4) Solving each subtask in parallel using the gradient descent method; 5) Determining the final identification result by combining the objective function value, the sparsity of the identified load, and the number of damaged elements. This invention is applicable to parallel computation, enabling simultaneous identification of structural dynamic loads and local damage, and is suitable for solving the overall sparsity problem of simultaneous identification of dynamic loads and local damage over long time spans.

Description

Technical Field

[0001] This invention belongs to the field of parallel computing technology, and in particular relates to a method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel. Background Technology

[0002] For engineering structures subjected to dynamic loads, dynamic loads are the primary external factor influencing the structure's dynamic behavior. Accurate monitoring of structural dynamic loads is beneficial for operation and maintenance management and is crucial for ensuring the safe operation of the structure. In practice, commonly used dynamic load acquisition methods can be divided into two main categories: direct measurement using force sensors and indirect inversion based on structural response. Direct identification methods typically involve installing force sensors at the location of the dynamic load, directly acquiring the structural dynamic load without complex algorithms. However, due to factors such as the operating environment, equipment costs, and installation difficulties, direct measurement methods are sometimes difficult to apply in practice. Unlike direct measurement methods, indirect measurement methods utilize measured structural responses, such as acceleration, displacement, and strain, combined with the structure's dynamic mapping characteristics to construct load identification equations, which are then solved using specific algorithms. Indirect measurement methods have unique advantages over direct measurement methods and can be widely applied to various conditions where dynamic loads are difficult to measure directly, such as wind loads and bridge moving loads.

[0003] Besides structural dynamic loads, changes in the dynamic parameters of the engineering structure itself are also a significant factor affecting its dynamic behavior. Since entering service, engineering structures are prone to damage and accumulation due to adverse factors such as environmental erosion, material aging, and long-term dynamic effects, which alters their dynamic characteristics. To promptly detect structural damage and ensure operational safety, damage identification research has received considerable attention. Past research has proposed numerous structural damage identification methods to address different practical needs, such as damage location and quantitative identification methods based on dynamic fingerprint changes and data-driven damage location identification methods. Among existing damage identification strategies, those based on simultaneously inverting external loads and structural damage based on structural response are of paramount importance. This is because, in practice, structural dynamic loads and structural damage are often unknown and difficult to directly perceive and detect.

[0004] In the research and application of simultaneous identification of structural dynamic loads and damage, the simultaneous identification strategy based on model correction has attracted much attention. This type of method considers the simultaneous identification problem as a parameter optimization problem, using the difference between the model-calculated response and the measured response of the structure, and constructing an optimization objective function through a certain mathematical model (such as the 2-norm) and solving it. Here, the optimization variables are usually selected as the parameters to be identified, namely the structural dynamic load and structural damage. For example, Zhu et al. [Zhu H, Mao L, Weng SA sensitivity-based structural damage identification method with unknown input excitation using transmissibility concept[J]. Journal of Sound and Vibration,2014,333(26):7135-7150] proposed a method for simultaneous identification of dynamic loads and damage based on damage sensitivity analysis. This method groups the measurement points and uses the load regularization identification method to construct a linear relationship between the two groups of measurement points. Then, by iteratively correcting the damage factor of each unit of the structure, the simultaneous identification of load and damage can be achieved. Meanwhile, a Chinese invention patent (authorization announcement number: CN110017929B) discloses a method for simultaneous identification of ship-collision load and damage to bridges based on substructure sensitivity analysis. One of the key steps in this strategy is to simultaneously identify ship-collision load and damage using a model correction method based on sensitivity analysis. Furthermore, a Chinese invention patent (authorization announcement number: CN106202789B) discloses a collaborative identification method for bridge moving load and damage. This method introduces sparse regularization based on model correction to ensure that bridge damage only occurs in local areas of the structure.

[0005] Currently, in research and application, methods for simultaneously identifying structural dynamic loads and structural damage, developed based on model modification, have achieved good results. However, because simultaneous identification of loads and damage requires identifying many parameters, and the iterative process requires updating the structural finite element model and calculating the corresponding system mapping matrix, the aforementioned methods are prone to drawbacks such as long computation time and easy misjudgment of damage when facing conditions with large load time spans and many discrete structural elements. Summary of the Invention

[0006] The purpose of this invention is to provide a method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel. By pre-dividing structural units into a finite number of groups, the original problem of simultaneous identification of loads and damages is decomposed into a finite number of independently runnable problems with fewer optimization parameters for simultaneous identification of dynamic loads and local damages. Furthermore, it integrates two acceleration strategies: dynamic load regularization identification and parallel computing, to achieve a low-time solution for simultaneous identification of structural dynamic loads and damages. On the one hand, it ensures the sparsity of the damage identification results in physical space, and on the other hand, it ensures that the aspects disclosed in this invention are applicable to parallel computing.

[0007] To achieve the above objectives, the present invention provides the following technical solution:

[0008] A method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel includes the following steps:

[0009] S1. Accelerometers are placed on the surface of the beam structure to collect the structural response, and the structural response is stored in a two-dimensional matrix form in combination with the moving time window.

[0010] S2. Use beam elements to build a finite element model, and divide the beam elements into n... eg Group n, construct n eg Each identification subtask can operate independently;

[0011] S3. For each identification subtask, the beam element damage factor contained in the group is used as the optimization parameter, and the optimization objective function is constructed by combining matrix regularization.

[0012] S4. Use gradient descent to solve each recognition subtask in parallel;

[0013] S5. Comprehensive objective function, identification of load sparsity, selection of dynamic load and structural damage final identification results.

[0014] The specific steps for arranging acceleration sensors on the surface of the beam structure to collect the structural response, and storing the structural response in a two-dimensional matrix form using a moving time window, are as follows:

[0015] S101. Acceleration sensors are uniformly arranged on the surface of the beam structure, and the acceleration sensors are used to record the acceleration time history response of the structure.

[0016] S102. Combining the moving time window, the acceleration time history response of the recorded structure is expressed in matrix form as follows:

[0017]

[0018] Where B represents a matrix containing the response information of all measurement points, and the subscript N s This indicates that the number of measuring points is N. sOne, B i Represented as the response to the i-th measurement point, the B i Represented as:

[0019]

[0020] Where w (italicized) represents the number of time windows, k (italicized) is the number of sampling points in each time window, k0 is the number of sampling points in the overlapping part of two adjacent time windows, and b i This represents the acceleration response corresponding to the i-th (italicized) sampling point.

[0021] Furthermore, the finite element model is established using beam elements, and the beam elements are divided into n... eg Group n, construct n eg Each identification subtask can be run independently, specifically:

[0022] S201. Use 2-node, 4-DOF beam elements to establish a finite element model of the beam structure;

[0023] S202, Divide the beam elements of the structure into n eg The number of beam elements contained in each group is adjustable, and the adjustableness allows a beam element to be included in different groups;

[0024] S203. For each beam element, a subtask for simultaneous load and damage identification is constructed in groups, requiring the generation of n subtasks in total. eg Each identification subtask can be computed in parallel. Each identification subtask only considers the beam elements contained in the group to be damaged, while the remaining beam elements are healthy.

[0025] Furthermore, for each identification subtask, the optimization objective function is constructed using the unit damage factor contained in the group as the optimization parameter, combined with matrix regularization, as follows:

[0026] S301. Structural damage is simulated using element stiffness reduction, and a dictionary matrix is ​​introduced to represent the dynamic loads. The first m (italicized) modal approximations of the structure's dynamic behavior are used to describe the structure's dynamic behavior. Modal displacements and modal velocities are used to represent the structure's vibration state at time 0. Specifically:

[0027]

[0028] Among them, K d K represents the stiffness matrix after structural damage. ei Let θ represent the stiffness matrix of the i-th element, N represent the number of structural elements, and θ represent the stiffness matrix of the i-th element. i ∈[0,1) represents the damage factor of the i-th unit, reflecting the degree of damage to the unit. θ is used for identification. i The range can be narrowed down, such as θ i∈[0,0.9];

[0029] Secondly, a moving time window is used to divide the time domain of the problem to be identified into n. tw For each short time window, an identification equation is constructed. When analyzing any given time window, the start time is set to 0, and the window length is T. w The time domain corresponding to the identification problem can be denoted as (0, T). w The j-th load within a time window is expanded using a proof (matrix) dictionary, denoted as:

[0030]

[0031] Where, vector f j Let D represent the discrete form of the j-th load, where D = [d1, d2, ..., dj]. n ] represents a matrix dictionary, whose atoms d i Satisfy ||d i || = 1, index n t This indicates the number of atoms contained in the load dictionary (matrix dictionary), α j Let represent the participation coefficients corresponding to the atoms in the matrix dictionary used for the j-th load; using the first m modal approximations to describe the dynamic behavior of the structure, the initial vibration state of the structure at time 0 can be described as:

[0032]

[0033] In the formula y Nm and These represent the initial modal displacement and initial modal velocity of the structure at order m, respectively.

[0034] Assuming the current subtask corresponds to a unit group containing n units, and their unit numbers are consecutive, then the current structural damage state can be denoted as θ = [0,…,0,θ]. i ,θ i+1 …,θ i+n-1 ,0,…,0] T When grouping units, it is important to ensure that i ≥ 1 and (i + n - 1) ≤ N.

[0035] S302. Based on the time-invariant and linear superposition principles of dynamic load application location, the relationship between structural dynamic load, initial conditions, and structural response is organized and expressed as follows:

[0036]

[0037] Among them, B i H represents the matrix form of the response at the i-th measurement point of the structure. ij(θ) represents the structural mapping matrix between the j-th load and the i-th measurement point, calculated by the modal superposition method, H iY (θ) represents the mapping matrix between the initial vibration of the structure and the i-th measuring point, which can also be calculated by the modal superposition method. The subscript N s and N f These represent the number of measuring points and the number of loads, respectively. Y represents the participation coefficient corresponding to the atom in the dictionary used by the i-th load within the j-th time window. j H represents the initial conditions corresponding to the j-th time window, and the mapping matrix is ​​H. f (θ) represents the part corresponding to the load, H y (θ) represents the part related to the initial vibration of the structure, F represents the dynamic load information of the structure, Y represents the initial condition information of the structure, and H f (θ) and H y In (θ), some elements may have large numerical differences. Before establishing the optimized recognition equation, the mapping matrix H is first... y (θ) is adjusted, and the adjusted mapping matrix is ​​denoted as H. y adjust (θ), any column of which, for example, the i-th column, is calculated by the following formula:

[0038]

[0039] In the formula, [:,i] represents taking the i-th column of the matrix, and C represents a positive constant, defined as Here, Q represents matrix H. f The number of columns in (0); using the adjusted mapping matrix, the response of the structure under the combined action of dynamic load and initial conditions can be further rewritten by formula (6) as:

[0040]

[0041] Where P = [F, Y] adjust ] represents the structural excitation matrix, F represents the structural dynamic load information, Y adjust This represents the method for representing the initial vibration state of the structure after coefficient adjustment, where A(θ) represents the system mapping matrix and θ represents the structural damage state.

[0042] When the damage state of the structure is θ, the corresponding dynamic load and initial conditions are calculated by the following formula:

[0043]

[0044] Where P(i,j) represents the element corresponding to the i-th row and j-th column of matrix P, F e (θ) and Y e(θ) represent the load identification value and the initial condition identification value, respectively, and λ represents the regularization parameter;

[0045] S303, the problem of simultaneously identifying dynamic loads and damages in optimized structures, is represented as follows:

[0046]

[0047] The objective function g(θ) is defined as follows:

[0048]

[0049] Furthermore, the partial derivative calculations involved in the gradient descent method are approximated using the difference method, and S4 specifically involves:

[0050] S401. Set initial values ​​for the algorithm parameters, including the regularization parameter λ, used to calculate the mapping matrix H. y adjust The adjustment coefficient C of (θ) is set, the damage vector update step size step and the loop variable itt = 2;

[0051] S402. Set the initial value of the damage vector as: θ(1)=[0,…,0,θ i (1),θ i+1 (1)…,θ i+n-1 (1),0,…,0] T Combined with formula (10), the objective function g[θ(1)] is calculated;

[0052] S403. For the itt-th iteration and itt≥2, calculate the gradient vector J of the objective function, specifically:

[0053]

[0054] The calculation of partial derivatives involved in the calculation process can be approximated by the finite difference method;

[0055] S404. Calculate the updated value of the damage vector, taking the temporary variable θ. temp Specifically:

[0056]

[0057] For θ temp Apply constraints such that if any element is greater than 0.9, then that element is set to 0.9; if any element is less than 0, then that element is set to 0.

[0058] S405. Determine whether the updated variable is acceptable. If g[θ] temp≥g[θ(itt - 1)], then this temporary update is not accepted, take θ(itt) = θ(itt - 1), g[θ(itt)] = g[θ(itt - 1)], reduce the step size by setting step = 0.5×step, increment the loop variable and set itt = itt + 1; if g temp < g[θ(itt - 1)], then this temporary update is accepted, take θ(itt) = θ temp , g[θ(itt)] = g[θ temp , enlarge the step size by setting step = 2×step, increment the loop variable and set itt = itt + 1;

[0059] S406. Determine whether the iteration meets the stopping condition. If it does, end the solution; if not, determine whether the damage variable has been updated in this iteration step. If there is an update, jump to S403; if not, jump to S404.

[0060] Further, updating the damage factor of the beam element only updates the beam elements included in the beam element group where the current subtask is located, and the remaining beam elements not included are defaulted to be healthy.

[0061] Further, the comprehensive objective function value represents the square of the Frobenius norm of the difference between the structural calculation response and the measured response, the identification load sparsity represents the number of non - zero elements in the identified structural excitation matrix, and the number of damaged elements represents the number of finite - element elements with a damage factor greater than 0 in the damage identification result.

[0062] Further, the specific steps of S5 are as follows:

[0063] S501. Compare the final objective function values of each identification subtask. The maximum and minimum values among all the final objective function values are g max and g min , and the difference between them is Δg;

[0064] S5, cyclic through the subtasks. The subtask cycle is that if the final objective function of a subtask is less than or equal to (g min + εΔg), then this subtask is considered to have a smaller objective function and can be selected for the next judgment, where ε is a fixed constant greater than 0 and less than 1; conversely, if the final objective function of a subtask is greater than (g min + εΔg), then the identification result of this subtask is directly discarded;

[0065] S503. For all subtasks with an objective function less than or equal to (g min + εΔg), calculate the judgment parameter corresponding to each task in combination with the Bayesian information criterion. The formula is:

[0066]

[0067] Where bic represents the judgment parameter calculated based on the Bayesian information criterion, n represents the number of elements in the response matrix B, and k nnz The activation matrix P represents the recognition. e The number of non-zero elements and the damage identification result θ e The sum of the number of non-zero elements in the middle;

[0068] S504, Compare the objective function to be less than or equal to (g min For all subtasks of +εΔg), the damage identification result corresponding to the subtask with the smallest bic value is taken as the final damage identification result;

[0069] S505. Calculate the structural dynamic load based on the final damage identification result of S504. If a certain moment is included in multiple time windows, the load identification result at that moment is taken as the average value of the identification results corresponding to multiple time windows.

[0070] Furthermore, S3 and S4 are parallelizable components.

[0071] Furthermore, parallel computation can be achieved by simultaneously using S3 and S4 on a parallel platform to solve multiple recognition subtasks.

[0072] The beneficial technical effects of the present invention are at least as follows:

[0073] The present invention provides a method for simultaneous identification of structural loads and local damage that can be computed in parallel. It integrates two acceleration strategies: dynamic load regularization identification and parallel computing. The method discloses a method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel. This method realizes the parallel solution of simultaneous identification of dynamic loads and local damage, and is suitable for low-time sparse solution of large-scale dynamic load and damage synchronous identification problems. Attached Figure Description

[0074] The present invention will be further described with reference to the accompanying drawings, but the embodiments in the drawings do not constitute any limitation on the present invention. For those skilled in the art, other drawings can be obtained based on the following drawings without creative effort.

[0075] Figure 1 This is a schematic diagram of the main implementation process of the method for simultaneous identification of structural loads and local damage that can be calculated in parallel according to an embodiment of the present invention.

[0076] Figure 2 This is a schematic diagram of the parallel analysis mode in the method for simultaneous identification of structural loads and local damage that can be calculated in parallel according to an embodiment of the present invention.

[0077] Figure 3This is a schematic diagram of the cantilever beam experimental model used in an embodiment of the present invention.

[0078] Figure 4 This diagram illustrates the acceleration response at the free end of a cantilever beam.

[0079] Figure 5 This diagram illustrates the comparison between simulated loads and identified loads.

[0080] Figure 6 This diagram illustrates the comparison between simulated damage and identified damage.

[0081] Figure labels: 1-Acceleration sensor; 2-Structural dynamic load. Detailed Implementation

[0082] Embodiments of the present invention are described in detail below. Examples of these embodiments are shown in the accompanying drawings, wherein the same or similar reference numerals denote the same or similar elements or elements having the same or similar functions throughout. The embodiments described below with reference to the accompanying drawings are exemplary and are only used to explain the present invention, and should not be construed as limiting the present invention. Specific implementation method one:

[0084] like Figure 1-2 As shown in the figure, the method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel according to the embodiments of the present invention includes the following main steps:

[0085] Step 1: Arrange accelerometers on the surface of the beam structure to collect the structural response. Combined with a moving time window, store the structural response in a two-dimensional matrix. For the i-th measurement point, use a time window to capture the structural response. The data captured by each time window are stored in a matrix format.

[0086]

[0087] Where w represents the number of time windows, k is the number of sampling points in each time window, k0 is the number of sampling points in the overlapping part of two adjacent time windows, and b i The acceleration response corresponding to the i-th sampling point is represented by the response matrix formed by multiple sampling points, which is:

[0088]

[0089] Where B represents a matrix containing the response information of all measurement points, and the subscript N s This indicates that the number of measuring points is N. s One, B i This is represented as the response to the i-th measurement point.

[0090] Step 2: Establish a finite element model using beam elements, and divide the beam elements into n... eg Group n, construct neg Each identification subtask can operate independently.

[0091] First, a finite element model of the beam structure is established using 2-node, 4-DOF beam elements; second, the beam elements of the structure are divided into n... eg There are n groups, each containing an adjustable number of elements, allowing a single element to be included in different groups; then, each element group constructs a load and damage simultaneous identification subtask, generating a total of n... eg Each subtask is a parallelizable identification task; here, each subtask considers only the units contained in its group to be potentially damaged, while the remaining units are healthy; the above subtasks are mutually independent.

[0092] Step 3: For each subtask, construct the optimization objective function by using the unit damage factor contained in the group as the optimization parameter and combining it with matrix regularization.

[0093] The analysis process for each subtask is as follows:

[0094] First, structural damage is simulated using element stiffness reduction, specifically:

[0095]

[0096] Among them, K d K represents the stiffness matrix after structural damage. ei Let θ represent the stiffness matrix of the i-th element, N represent the number of structural elements, and θ represent the stiffness matrix of the i-th element. i ∈[0,1) represents the damage factor of the i-th unit, reflecting the degree of damage to the unit. θ is used for identification. i The range can be narrowed down, such as θ i ∈[0,0.9]; secondly, a moving time window is used to divide the time domain of the problem to be identified into n. tw For regions with relatively short durations, overlap between adjacent time windows is allowed. An identification equation is constructed for each time window. Without loss of generality, when analyzing any given time window, the starting time of the currently considered time window can be set to 0, and the window length to T. w Then the time domain corresponding to the identification problem can be denoted as (0, T). w The j-th load within a time window is represented by a dictionary expansion:

[0097]

[0098] Where, vector f j Let D represent the discrete form of the j-th load, where D = [d1, d2, ..., dj]. n ] represents the load dictionary, whose atoms d i Satisfy ||d i || = 1, index nt Indicates the number of atoms contained in the load dictionary used, α j Let represent the participation coefficients corresponding to the atoms in the dictionary used for the j-th load; using the first m modal approximations to describe the dynamic behavior of the structure, the initial vibration state of the structure at time 0 can be described as:

[0099]

[0100] Among them, y Nm and Let θ represent the initial modal displacement and initial modal velocity of the structure at order m, respectively. Assuming that the element group corresponding to the current subtask contains n elements and their element numbers are consecutive, the current structural damage state can be denoted as θ.

[0101] =[0,…,0,θ i ,θ i+1 …,θ i+n -1,0,…,0] T When grouping elements, it is important to ensure that i ≥ 1 and (i + n - 1) ≤ N. Considering that the location of the dynamic load remains unchanged, based on the principle of linear superposition, the structural response under the combined action of the dynamic load and initial conditions can be organized as follows:

[0102]

[0103] Among them, B i H represents the matrix form of the response at the i-th measurement point of the structure. ij (θ) represents the structural mapping matrix between the j-th load and the i-th measuring point, which can be calculated using the modal superposition method. H iY (θ) represents the mapping matrix between the initial vibration of the structure and the i-th measuring point, which can also be calculated by the modal superposition method; subscript N s and N f These represent the number of measuring points and the number of loads, respectively. Y represents the participation coefficient corresponding to the atom in the dictionary used by the i-th load within the j-th time window. j H represents the initial conditions corresponding to the j-th time window, and the mapping matrix is ​​H. f (θ) represents the part corresponding to the load, H y (θ) represents the part related to the initial vibration of the structure, F represents the dynamic load information of the structure, Y represents the initial condition information of the structure, and H f (θ) and H y In (θ), some elements may have large numerical differences. Before establishing the optimized recognition equation, the mapping matrix H is first... y (θ) is adjusted, and the adjusted mapping matrix is ​​denoted as H. y adjust(θ), any column of which, for example, the i-th column, is calculated by the following formula:

[0104]

[0105] In the formula, [:,i] represents taking the i-th column of the matrix, and C represents a positive constant, defined as Here, Q represents matrix H. f The number of columns in (0); using the adjusted mapping matrix, the response of the structure under the combined action of dynamic load and initial conditions can be further rewritten by formula (6) as:

[0106]

[0107] Where P = [F, Y] adjust ] represents the structural excitation matrix, Y adjust Represents the mapping matrix H y adjust The method for representing the initial vibration state of the structure corresponding to (θ), the system mapping matrix A(θ) = [H f (θ),H y adjust Finally, based on formula (8) and combined with the matrix regularization identification method of dynamic load, it can be seen that when the damage state of the structure is θ, the corresponding dynamic load and initial conditions can be calculated by the following formula:

[0108]

[0109] Where P(i,j) represents the element corresponding to the i-th row and j-th column of matrix P, F e (θ) and Y e (θ) represent the load identification value and the initial condition identification value, respectively, and λ represents the regularization parameter, which can be appropriately selected according to the Bayesian Information Criterion (BIC). Using formula (9), the objective function g(θ) is further defined as the square of the Frobenius norm of the difference between the calculated response and the measured response of the structure, mathematically denoted as:

[0110]

[0111] The problem of simultaneously identifying structural dynamic loads and damages can be reduced to solving the following optimization problem:

[0112]

[0113] Where, θ e Values ​​indicating structural damage;

[0114] Step 4. Use the gradient descent method to solve each sub-task, that is, solve formula (11). The specific process is as follows:

[0115] Start: Step (A). Set the initial values of the algorithm parameters, including the regularization parameter λ, the adjustment coefficient C for calculating the mapping matrix H yadjust (θ); the update step size step of the damage vector; set the loop variable itt = 2;

[0116] Step (B). Set the initial value of the damage vector as: θ(1) = [0, …, 0, θ i (1), θ i+1 (1) …, θ i+n-1 (1), 0, …, 0] T , and calculate the objective function g[θ(1)] in combination with formula (10);

[0117] Step (C). For the itt (itt ≥ 2) -th iteration, calculate the gradient vector J of the objective function, specifically:

[0118]

[0119] The partial derivative calculation involved in the calculation process can be approximately replaced by the difference method;

[0120] Step (D). Calculate the updated value of the damage vector, and take the temporary variable θ temp as:

[0121]

[0122] Constrain θ temp . If a certain element is greater than 0.9, then this element is taken as 0.9; if a certain element is less than 0, then this element is taken as 0;

[0123] Step (E). Judge whether the updated variable is accepted; if g[θ temp ≥ g[θ(itt - 1)], then this temporary update is not accepted, take θ(itt) = θ(itt - 1), g[θ(itt)] = g[θ(itt - 1)], reduce the step size and set step = 0.5 × step, increase the loop variable and take itt = itt + 1; if g[θ temp < g[θ(itt - 1)], then this temporary update is accepted, take θ(itt) = θ temp , g[θ(itt)] = g[θ temp , enlarge the step size and set step = 2 × step, increase the loop variable and take itt = itt + 1;

[0124] In step (F), determine whether the iteration meets the stopping condition. If it does, end the solution; if it does not, determine whether the damage variable has been updated in this iteration step. If it has been updated, jump to step (C); if it has not been updated, jump to step (D).

[0125] Note that during the above solution process, the damage factor that needs to be updated is only the unit included in the unit group used by the current subtask; other units not included are considered healthy by default.

[0126] Step 5: Combine the objective function value, identified load sparsity, and number of damaged elements, and use the Bayesian information criterion to determine the final identification result. Here, the objective function value is represented by the square of the Frobenius norm of the difference between the calculated and measured structural response; the identified load sparsity is represented by the number of non-zero elements in the identified structural excitation matrix; and the number of damaged elements represents the number of finite element elements with a damage factor greater than 0 in the damage identification result. Step 1: Compare the final objective function values ​​of each subtask, and denote the maximum and minimum values ​​among all final objective function values ​​as gmax and gmin, respectively. max and g min Let the difference between the two be Δg; for the subtask loop, if the final objective function of a subtask is less than or equal to (g) min If the objective function of a subtask is greater than or equal to εΔg, then the subtask is considered to have a smaller objective function and can proceed to the next step of judgment. Here, ε is a fixed constant greater than 0 and less than 1, representing the acceptable fluctuation range of the objective function. Conversely, if the final objective function of a subtask is greater than or equal to εΔg, then the subtask is considered to have a smaller objective function and can proceed to the next step of judgment. min If the objective function is greater than or equal to (g), then the subtask is considered to have a larger objective function, and the result of the subtask is directly discarded. The second step is to handle subtasks with an objective function less than or equal to (g). min For all subtasks of +εΔg), the decision parameters for each task are calculated using the Bayesian Information Criterion (BIC), as shown in the following formula:

[0127]

[0128] Where bic represents the judgment parameter calculated based on the Bayesian information criterion, n represents the number of elements in the response matrix B, and k nnz The activation matrix P represents the recognition. e The number of non-zero elements and the damage identification result θ e The sum of the number of non-zero elements; compared to the objective function being less than or equal to (g min For all subtasks of +εΔg), the damage identification result corresponding to the subtask with the smallest bic value is taken as the final damage identification result; further, the identification result of the structural dynamic load is calculated by formula (9) and formula (4). If a certain moment is included in multiple time windows, the load identification result at that moment is taken as the average value of the identification results corresponding to multiple time windows. Specific Implementation Method Two:

[0130] As attached Figure 3 As shown, this embodiment is based on Specific Implementation Method 1. Consider a cantilever beam model under load. The cantilever beam is 0.6m long and is uniformly divided into 12 beam elements, each containing 2 nodes and 4 degrees of freedom. The bending stiffness of the model is EI = 15.96 Nm. -2 Linear density ρA = 1.062 kg m -1 Six accelerometers were placed at distances of 0.1m, 0.2m, 0.3m, 0.4m, 0.5m, and 0.6m from the support. The influence of the added mass of the sensors was considered during finite element modeling, with each sensor weighing 0.039kg. The first three natural frequencies of the undamaged model were 4.92Hz, 30.72Hz, and 85.73Hz. Now, a 15% stiffness reduction is assumed in element 3, and a load is applied at a distance of 0.35m from the support. The time history signal of the load is as follows:

[0131]

[0132] Using the first eight modal data of the structure, the modal superposition method was selected to simulate the structural vibration response. The damping ratios of the first eight modes were set to 2.05%, 1.2%, 1.12%, 0.82%, 0.77%, 0.5%, 0.77%, and 0.5 / 100, respectively. The time discretization interval for the forward analysis was 1 / 5000 s. Six accelerometers were uniformly arranged on the lower surface of the cantilever beam, with a sampling frequency of 1000 Hz. The acceleration response was simulated using white noise, taking into account the influence of measurement noise, with a signal-to-noise ratio of 30 dB for each measurement point. The sampling time at each measurement point was 4 seconds. The acceleration response at a distance of 0.6 m is shown in the appendix. Figure 4 .

[0133] Based on the simulated measured responses at six acceleration measurement points, the specific implementation steps for simultaneous identification of dynamic load and local damage disclosed in this invention are as follows:

[0134] (1) An accelerometer was used to collect the structural response of the cantilever beam at distances of 0.1m, 0.2m, 0.3m, 0.4m, 0.5m, and 0.6m from the support. The sampling frequency was 1000Hz, and the sampling time was 4s. Time windows were selected as (0s, 0.25s], (0.25s, 0.5s], (0.5s, 0.75s], (0.75s, 1s], (1s, 1.25s], and (1.25s, 1... The response of the measurement point was captured in 16 time windows: (1.5s, 1.75s), (1.75s, 2s), (2s, 2.25s), (2.25s, 2.5s), (2.5s, 2.75s), (2.75s, 3s), (3s, 3.25s), (3.25s, 3.5s), (3.5s, 3.75s), and (3.75s, 4s). The response was then stored in columns as a response matrix B.

[0135] (2) Establish the structural finite element model, as shown in the attached figure. Figure 3 As shown; here, the unit numbering starts from the support end and is numbered along the positive x-axis as Unit 1, Unit 2, Unit 3, Unit 4, Unit 5, Unit 6, Unit 7, Unit 8, Unit 9, Unit 10, Unit 11, Unit 12 respectively; based on the unit numbering, this embodiment will be divided into 23 groups, namely: [Unit 1], [Unit 2], [Unit 3], [Unit 4], [Unit 5], [Unit 6], [Unit 7], [Unit 8], [Unit 9], [Unit 1], [Unit 12], [Unit 13], [Unit 14], [Unit 15], [Unit 16], [Unit 17], [Unit 18], [Unit 19], [Unit 10], [Unit 11], [Unit 12 ...1], [Unit 12], [Unit 13], [Unit 14], [Unit 15], [Unit 16], [Unit 17], [Unit 18], [Unit 19], [Unit 10], [Unit 11], [Unit 12], [Unit 13], [Unit 14], [Unit 15], [Unit 16], [Unit 17], [Unit 18], [Unit 19], [Unit 10], [Unit 11], [Unit 12], [Unit 13], [Unit 14], [Unit 15], [Unit 16], [Unit 17], [Unit 18], [Unit 1 [Unit 10], [Unit 11], [Unit 12], [Unit 1, Unit 2], [Unit 2, Unit 3], [Unit 3, Unit 4], [Unit 4, Unit 5], [Unit 5, Unit 6], [Unit 6, Unit 7], [Unit 7, Unit 8], [Unit 8, Unit 9], [Unit 9, Unit 10], [Unit 10, Unit 11], [Unit 11, Unit 12]; Based on the 23 groups, construct 23 independently computable parallel subtasks.

[0136] (3) Load dictionary is constructed using load shape functions. The idea behind constructing load shape functions is to treat the time history curve of the load in the time domain as the deflection of a "time beam". In this embodiment, the 0.25s duration is divided into 50 "time beam elements" using 51 nodes, with each "time beam element" having a duration of 0.005s. The load time history curve is fitted using the shape functions of a 2-node, 4-DOF beam element. A total of 51 × 2 = 102 load shape functions are required for the 0.25s time history. Assuming no structural damage occurs, the mapping matrix H is calculated. y adjustThe adjustment coefficient C of (θ) is used to generate A(0). When calculating A(0), the first 8 modal information of the structure is selected, and the discrete interval of the analysis time is 1 / 5000s. Secondly, based on A(0) and the response matrix B, the FISTA (FastIterative Shrinkage Thresholding Algorithm) algorithm is used to iteratively solve formula (9) to realize the load estimation. The regularization parameter λ is selected using the Bayesian information criterion. Finally, keeping the constant C and the regularization parameter λ unchanged, 23 optimization equations for simultaneous identification of load and damage are constructed in parallel according to the 23 unit groups.

[0137] (4) Gradient descent is used to solve each subtask in parallel. The parallel computing platform consists of 8 computing nodes and is built using the Parallel Computing Toolbox in Matlab. Parallel computing is programmed using the parfor mode. For each subtask, the initial iteration value is set to structural health, and the iteration step size is set to step = 0.05. The partial derivatives are calculated using the difference approximation, as shown in the following formula:

[0138]

[0139] In the formula, Δθ=[0,…,0,θ i =δθ,0,…,0] T δθ is set to 0.005. The shutdown condition is set to step ≤ 0.00001 or the number of iterations reaches 50.

[0140] (5) Compare the objective functions after iterative solutions of the 23 sub-tasks, and eliminate the sub-task identification results with larger objective functions by setting ε to 5%. For the selected sub-task identification results, use the Bayesian criterion to select the sub-task identification result with the smallest bic value as the final result, and calculate the structural dynamic load identification result by combining formula (9) and formula (4).

[0141] Appendix Figure 5 The simulation of dynamic loads and the identification of dynamic loads were compared. Figure 6 The simulated structural damage from the examples and the response-based damage identification results were compared. (From the appendix...) Figure 5 and attached Figure 6As can be seen, a parallel computational method for simultaneous identification of structural loads and local damage can effectively identify simulated dynamic loads and structural damage with high accuracy. It is noted that in this embodiment, during damage identification, the elements are grouped in physical space, ensuring significant sparsity in the damage identification results and significantly reducing the variables requiring optimization during the solution process, thus accelerating algorithm convergence. Furthermore, the identification subtasks formed by the element groupings are mutually independent, making this invention suitable for parallel computation and reducing computational time. Finally, when establishing the system matrix, this invention only needs to consider a window length of 0.25s, much smaller than the total analysis time of 4s. Therefore, this invention can handle the problem of simultaneous identification of dynamic loads and local damage over a long time span using a smaller-dimensional system matrix.

[0142] Although embodiments of the invention have been shown and described, those skilled in the art will understand that various changes, modifications, substitutions and variations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the claims and their equivalents.

Claims

1. A method for simultaneous identification of structural dynamic loads and local damages which can be calculated in parallel, characterized in that, It includes the following steps: S1. Accelerometers are placed on the surface of the beam structure to collect the structural response, and the structural response is stored in a two-dimensional matrix form in combination with the moving time window. S2, a finite element model is established using beam elements, and the beam elements are divided into n eg n eg independent recognition sub-tasks are constructed​ S3. For each identification subtask, the beam element damage factor contained in the group is used as the optimization parameter, and the optimization objective function is constructed by combining matrix regularization. Specifically, for each identification subtask, the optimization objective function is constructed using the unit damage factor contained in the group as the optimization parameter and matrix regularization. The specific steps are as follows: S301, simulate structure damage using unit stiffness reduction and introduce dictionary matrix to represent structure dynamic load, use structure pre m order modal to describe structure dynamic behavior, use modal displacement and modal velocity to represent structure vibration state at 0 time; S302. Based on the time-invariant and linear superposition principles of dynamic load application location, the relationship between structural dynamic load, initial conditions, and structural response is organized and expressed as follows: ； Where P=[F, Y] adjust ] represents the structural excitation matrix, F represents the structural dynamic load information, Y adjust This represents the method for representing the initial vibration state of the structure after coefficient adjustment, where A(θ) represents the system mapping matrix and θ represents the structural damage state. When the damage state of the structure is θ, the corresponding dynamic load and initial conditions are calculated by the following formula: ； Wherein, P( i , j ) represents the first digit of matrix P. i Line 1 j The element corresponding to the column, F e (θ) and Y e (θ) represent the load identification value and the initial condition identification value, respectively, and λ represents the regularization parameter; S303, the problem of simultaneously identifying dynamic loads and damages in optimized structures, is represented as follows: ； where the objective function g (θ) is defined as: ； S4. Use gradient descent to solve each recognition subtask in parallel; S5. Comprehensive objective function, identification of load sparsity, selection of dynamic load and structural damage final identification results.

2. The method of claim 1, wherein, The specific steps for arranging acceleration sensors on the surface of the beam structure to collect the structural response, and storing the structural response in a two-dimensional matrix form using a moving time window, are as follows: S101. Acceleration sensors are uniformly arranged on the surface of the beam structure, and the acceleration sensors are used to record the acceleration time history response of the structure. S102. Combining the moving time window, the acceleration time history response of the recorded structure is expressed in matrix form as follows: ； where B represents a matrix containing all the response information of the measuring points, the subscript N s denotes the number of measuring points N s is 1, B i represents the response to the i th measuring point, said B i is represented as: ； wherein, w denotes the number of time windows, k is the number of sampling points contained in each time window, k 0 is the number of sampling points of the overlapping part of two adjacent time windows, b i denotes the acceleration response corresponding to the i th sampling point.

3. The method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel according to claim 1, characterized in that, The beam unit is divided into n eg groups, and n eg independent identification sub-tasks are constructed, and the specific steps are as follows: S201. Establish a finite element model of the beam structure using 2-node, 4-DOF beam elements; S202, dividing the beam elements of the structure into n eg groups, the number of beam elements included in each group being adjustable, the adjustability allowing one beam element to be included in different groups; S203, each beam element group constructs a load and damage simultaneous identification subtask, a total of n eg parallel computing identification subtasks, each identification subtask only considers the damage of the beam elements contained in the group, and the rest of the beam elements are healthy.

4. The method of claim 1, wherein, The partial derivative calculations involved in the gradient descent method are approximated using the difference method. The specific steps in S4 are as follows: S401, set algorithm parameter initial value, including regularization parameter λ, used for calculating mapping matrix H y adjust adjustment coefficient of (θ) C， set damage vector update step size step and loop variable itt =2; S402、Set the initial value of the damage vector as: θ(1)=[0, …,0, θ i (1), θ i+1 (1)…, θ i+n-1 (1) , 0, …,0] T , combined with formula (10) to calculate the objective function g [θ(1)]; S403, Regarding the first itt Step iteration and itt ≥2, calculate the gradient vector J of the objective function; S404, calculate the update value of the damage vector, take the temporary variable θ temp ; S405. Determine whether the updated variable is acceptable; S406. Determine if the iteration meets the stopping condition. If it does, end the solution process. If the condition is not met, determine whether the damage variable has been updated in the iteration step. If it has been updated, jump to S403; otherwise, jump to S404.

5. The method of claim 1, wherein, Updating the beam element damage factor means only updating the beam elements included in the beam element group of the current subtask; other beam elements not included are considered healthy by default.

6. The method of claim 1, wherein, The comprehensive objective function value represents the square of the Frobenius norm of the difference between the calculated and measured structural responses; the identified load sparsity represents the number of non-zero elements in the identified structural excitation matrix; and the number of damaged elements represents the number of finite element elements with a damage factor greater than 0 in the damage identification results.

7. The method of claim 6, wherein, S5 specifically includes the following steps: S501. Compare the final objective function values ​​of each recognition sub-task. The maximum and minimum values ​​among all the final objective function values ​​are respectively... g max and g min The difference between the two is Δ g ; S502, a subtask loop is performed, and if a final target function of a subtask is less than or equal to ( g min + epsilon * delta g ), the subtask is considered to have a smaller target function, and a next step is selected to be entered, wherein epsilon is a fixed constant greater than 0 and less than 1; on the contrary, if the final target function of the subtask is greater than ( g min + epsilon * delta g ), the identification result of the subtask is directly discarded; S503、For all subtasks whose target functions are less than or equal to ( g min + ε Δ g ), the judgment parameters corresponding to each task are calculated in combination with the Bayesian information criterion, and the formula is as follows: ； wherein, bic represents a judgment parameter calculated based on Bayesian information criterion, n represents the number of elements in the response matrix B, k nnz represents the identified excitation matrix P e the number of non-zero elements in the damage identification result θ e the sum of the number of non-zero elements in the damage identification result θ S504, Compare the objective function to be less than or equal to ( g min +εΔ g Take all subtasks of ) bic The damage identification result corresponding to the subtask with the smallest value is taken as the final damage identification result; S505. Calculate the structural dynamic load based on the final damage identification result of S504. If a certain moment is included in multiple time windows, the load identification result at that moment is taken as the average value of the identification results corresponding to multiple time windows.

8. The method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel according to claim 1, characterized in that, Steps S3 and S4 are the parts that can be computed in parallel.

9. The method for simultaneous identification of structural dynamic loads and local damage that can be computed in parallel according to claim 8, characterized in that, Parallel computation means that multiple recognition subtasks can be solved simultaneously using steps S3 and S4 through a parallel platform.

Images