Bayesian optimization estimation method for kinetic parameters of rotating machinery

Abstract

The invention discloses a Bayesian optimization estimation method for kinetic parameters of a rotating machine, and the method comprises the steps: synchronously collecting a key phase signal and a transient vibration signal when the rotating speed of the rotating machine changes, and reconstructing a rotating speed pulse time sequence and a time domain signal; fourier transform is carried out on the vibration signal, and transient response data of a specified order is extracted; deleting the displacement offset in the transient response data; function fitting is carried out on the transient response data after offset elimination, and a transient response function with the rotating speed or the frequency as an independent variable and the amplitude as a dependent variable is generated; defining an optimization parameter space and a target function; and transmitting the target function and the optimization parameter space to a Bayesian optimization model to obtain an optimal parameter combination, and identifying the dynamic parameters of the rotating machinery based on the optimal parameter combination. The method can reduce the cost of computing resources, has wide parameter adaptability, and can estimate various kinetic parameters such as support stiffness and structural damping.

Description

Technical Field

[0001] The present application belongs to the technical field of rotating machinery dynamics, and specifically relates to a Bayesian optimization estimation method for rotating machinery dynamics parameters. Background Art

[0002] During the design and maintenance phases of rotating machinery, obtaining accurate rotor system dynamic parameters is crucial for precise dynamic simulation and optimal structural optimization. However, these dynamic parameters are influenced by numerous complex factors, including but not limited to structural layout, material selection, and installation method. Due to the diversity and complexity of these factors, accurate dynamic parameters cannot be directly obtained through theoretical calculations.

[0003] To model and simulate the dynamic characteristics of a rotor system, it is often necessary to combine limited measured data with known structural parameters and apply parameter estimation (identification) techniques to estimate the dynamic parameters of the rotor system. In engineering practice, modal testing is a commonly used method for obtaining the dynamic parameters of a rotor system. This technique obtains relevant modal data by conducting experiments while the mechanical system is at rest. This modal data is then processed and analyzed using specific formulas to calculate the dynamic parameters of the rotor system, such as support stiffness and structural damping.

[0004] However, modal testing technology has obvious limitations. On the one hand, since the test is performed when the mechanical system is stationary, it cannot reflect some of the dynamic characteristics of the rotating machinery during operation. For example, when the rotating machinery is running, dynamic factors such as bearing forces and oil film forces provided by the damper will have a significant impact on the dynamic characteristics of the rotor system. However, due to the limitations of the test conditions, modal testing technology cannot effectively capture these dynamic characteristics. On the other hand, modal testing technology has specific requirements for the test environment and requires specific boundary conditions, such as suspension and hoisting. In actual engineering applications, it is sometimes difficult to meet these specific test environment requirements. The above limitations make the current means of obtaining dynamic parameters insufficient in terms of accuracy and applicability, which in turn affects the simulation accuracy and optimization effect, and may lead to an extension of the R&D cycle and an increase in development costs for rotating machinery. Summary of the Invention

[0005] Given the shortcomings of existing modal testing techniques in terms of accuracy and applicability, this application aims to provide a Bayesian optimization method for estimating the dynamic parameters of rotating machinery. This method reduces computational resource overhead, exhibits broad parameter adaptability, and can estimate a variety of dynamic parameters, such as support stiffness and structural damping. Furthermore, when the derived parameters are used for rotor system dynamic simulation, the simulation model can simulate realistic vibration characteristics with guaranteed accuracy.

[0006] In order to achieve the above technical objectives, this application specifically adopts the following technical solutions:

[0007] In one aspect of the present application, a Bayesian optimization method for estimating dynamic parameters of a rotating machinery is provided, comprising the following steps:

[0008] S1: Synchronously collect the key phase signal S when the rotating machinery speed changes tacho And transient vibration signal x(t), and reconstruct the speed pulse time series and the time domain signal x clean (t);

[0009] S2: vibration signal x clean (t) performing Fourier transform to obtain transient spectrum data, and extracting transient response data of a specified order therefrom;

[0010] S3: deleting the displacement offset caused by shaft bending and / or installation deviation in the transient response data;

[0011] S4: Perform function fitting on the transient response data after eliminating the bias, and generate a transient response function g with the speed or frequency as the independent variable and the amplitude as the dependent variable. i (ω);

[0012] S5: Define the optimization parameter space and objective function:

[0013] The optimization parameter space includes at least one of support stiffness, structural damping parameter, rotor eccentricity, bearing force parameter and oil film force parameter;

[0014] Based on the transient response function g i (ω) Generate an objective function that receives the optimization parameters and calls the simulated rotor model to output a scalar value:

[0015]

[0016] Ω={ω sim (1),ω sim (2), …};

[0017] in, is the scalar value of the objective function f(x); N component ,N signal are the set number of frequency domain components and the set number of signals respectively; Ω represents the matrix |X sim All elements in the corresponding instantaneous speed vector, matrix |X sim| is the complex matrix of the specified order transient response extracted from the simulated vibration response data by Fourier transform in the vibration response analysis of the rotor model. The elements in the matrix represent the amplitude of the corresponding order vibration component at different frequencies; g i,j (Ω) represents the value of the transient response function of the i-th frequency domain component and the j-th signal at the speed Ω; |X sim | i,j (Ω) represents the i-th frequency domain component, the j-th signal corresponds to the instantaneous speed Ω |X sim |The value of the matrix; w i,j (Ω) represents the weight coefficient of the i-th frequency domain component and the j-th signal at the instantaneous speed Ω;

[0018] The smaller the scalar value, the closer the transient response obtained by the simulated rotor model is to the transient response function. When the scalar value shows a decreasing trend, it indicates that the degree of closeness between the transient response obtained by the simulated rotor model and the transient response function is higher.

[0019] S6: passing the objective function and the optimization parameter space to the Bayesian optimization model to obtain the optimal parameter combination, and identifying the dynamic parameters of the rotating machinery based on the optimal parameter combination.

[0020] In one embodiment, the denoised reconstructed tachometer pulse time series include:

[0021] S111: Get key phase signal S tacho The average value of the peak and valley values ​​T threshold as a trigger threshold;

[0022] S112: Get key phase signal S tacho Pulse front time series Calculate the pulse start time based on the trigger threshold

[0023]

[0024] Where Δt = 1 / f s , f s is the sampling frequency, i represents the index in the sequence;

[0025] S113: According to the predefined speed threshold ω threshold Delete the low-speed signals at the beginning and end of the signal to generate a valid data time series and speed sequence

[0026] S114: Speed ​​sequence Double smoothing optimization time series Generate a quadratically smoothed speed series

[0027] S115: Based on smooth speed sequence Correct the time series to generate the speed pulse time series

[0028] In one embodiment, the denoised reconstructed time domain signal x clean (t) include:

[0029] S121: Based on the speed pulse time series Constructing the corner area Use cubic spline interpolation to generate The corresponding angular domain time series

[0030] S122: Based on angular domain time series Mapping the transient vibration signal x(t) to the angular domain vibration signal x(θ);

[0031] S123: Use a sliding window to perform synchronous averaging on the angular domain vibration signal x(θ) to generate an angular domain signal x clean (θ);

[0032] S124: Based on the angle domain signal x clean (θ) generates the corresponding angular domain time series and time series in the time domain

[0033] S125: According to the frequency of the transient vibration signal x(t), the time series Uniform sampling generates time series t 5 ;

[0034] S126: The angle domain signal x clean (θ) is mapped to the time series t 5 Generate time domain signal x clean (t).

[0035] In one embodiment, step S2 includes:

[0036] S21: Based on the speed pulse time series and the angular domain vibration signal x(θ) to generate the angular domain signal x clean (θ), and according to the angular domain signal x clean (θ) corresponds to the time domain time series Construct the fundamental frequency cycle:

[0037]

[0038] Among them, N base is the set fundamental frequency denominator parameter;

[0039] S22: The time domain signal x clean (t) is mapped to the angular domain to generate the angular domain signal x clean (θ), and divided into blocks according to the fundamental frequency period;

[0040] S23: For each block of inner angle signal x clean (θ) performs a windowed fast Fourier transform and outputs a complex matrix X:

[0041] X=[X1,X2,…,X n ] T ,(n <f s / 2);

[0042] The i-th row in the matrix represents the rotation speed i / N base times the frequency domain component of the frequency vibration;

[0043] S24: extract the complex elements of the row of the specified order from the complex matrix X, calculate the amplitude and generate an amplitude sequence;

[0044] S25: Construct the angular velocity domain corresponding to the amplitude sequence, based on the time domain time series The instantaneous speed ω(i) is calculated by differential calculation and the instantaneous speed vector ω is generated by smoothing. 2 , forming the transient response data of amplitude-speed.

[0045] In one embodiment, step S3 includes: extracting amplitude values ​​from the transient spectrum data and storing them as a |X| matrix; for each row X in the |X| matrix, i Remove the offsets individually:

[0046]

[0047] Among them, X i,j is the element in the i-th row and j-th column of the |X| matrix; β is the attenuation factor; and T is the adaptive threshold.

[0048] In one embodiment, step S4 includes:

[0049] S41: Obtain the row vector of the amplitude sequence and the corresponding instantaneous speed vector ω based on the transient response data 2 ;

[0050] S42: With speed as the independent variable and amplitude as the dependent variable, a proxy function is constructed using cubic spline interpolation:

[0051]

[0052] S43: Minimize the following objective function and store the transient response function g i (ω):

[0053]

[0054] Where p represents the smoothing parameter; w j represents the weight function; represents the second-order derivative of the surrogate function; X i,j represents the element in the i-th row and j-th column of the amplitude sequence, i.e., the vibration amplitude corresponding to order i at the j-th instantaneous speed point; g i (ω) represents the predicted value of the surrogate function.

[0055] In one embodiment, defining the optimization parameter space in step S5 includes:

[0056] S51: defining optimization parameters and the value range of each optimization parameter;

[0057] S52: Convert each parameter to the normalized space through the function h(x):

[0058]

[0059] Where h(x) represents the logarithmic and normalized function; x represents the optimization parameter.

[0060] In one embodiment, defining the objective function in step S5 includes:

[0061] S53: performing inverse normalization and inverse logarithmization on the optimized parameters to restore the physical parameters;

[0062] S54: Input physical parameters into the rotor dynamics differential equation to update the rotor system mass matrix, stiffness matrix, damping matrix, gyro matrix, and transient matrix:

[0063]

[0064] Where y represents the node displacement of the rotor system; M, K, C, G, and N represent the mass matrix, stiffness matrix, damping matrix, gyro matrix, and transient matrix of the rotor system, respectively; and F represents the external force.

[0065] S55: Determine the vibration response y under the same measured speed / acceleration conditions sim , extract the response data matching the measurement point

[0066] S56: Response data Detrending is performed and Fourier transform is used to extract the transient response of the specified order to generate the simulation amplitude matrix |X sim |Matrix and speed vector ω sim ;

[0067] S57: Generate an objective function using a multi-peak Gaussian function weighting.

[0068] In one embodiment, the multi-peak Gaussian weight function is defined as:

[0069]

[0070] Among them, w(Ω) is the weight coefficient at the speed Ω, α i is the amplitude coefficient of the i-th Gaussian kernel, μ i is the central speed of the i-th Gaussian kernel, σ i is the standard deviation of the i-th Gaussian kernel used for weighted bandwidth control.

[0071] In one embodiment, step S6 includes:

[0072] S61: Read the optimization parameter space and calculate the objective function Establish a Gaussian process prior model:

[0073] f(x)~GP(m(x),k(x,x′;λ,l))

[0074] Where m(x) is the mean function; λ is the output scale hyperparameter that controls the output amplitude; l is the length scale hyperparameter that controls the smoothness of the function;

[0075] The covariance kernel function k(x,x') is calculated using the Matérn 5 / 2 kernel function:

[0076]

[0077] Where D represents the dataset;

[0078] S62: Based on the current observation data set Calculate the Gaussian process posterior distribution and select the next evaluation point x by the expected improvement acquisition function next :

[0079] α EI =E[max(f(x * )-f(x),0)]=∫max(f(x * )-f(x),0)p(f(x)|x)df(x)

[0080] Among them, E[·] is the expectation, x * is the optimal set of parameters currently observed, and p(·) represents the probability distribution;

[0081] S63: x next Substitute the objective function f(x) to calculate the scalar value L next , and the new data {x next ,L next}Update to dataset D to form a new dataset D';

[0082] S64: When the maximum number of evaluations or the convergence threshold is met, the current optimal parameter x is output. * .

[0083] The beneficial effects of this application are:

[0084] 1) This application will significantly reduce the requirements of the parameter estimation method on the test environment, has a wide range of parameter estimation, and can perform parameter estimation in high-dimensional parameter space.

[0085] 2) When applying this application for parameter estimation, engineers only need to collect vibration data from the rotating machinery rotor system during acceleration or deceleration, which reduces the requirements for the test environment. In addition, the method can apply vibration data from multiple points without significantly increasing computing resource overhead.

[0086] 3) This application can estimate a variety of parameters, including support stiffness, structural damping, bearing force parameters, damper oil film force parameters, etc. At the same time, changes in parameter types have no significant impact on the implementation of this method, that is, this method has high robustness to different parameters.

[0087] 4) This application can estimate multiple parameters at the same time and consider the coupling effect of each parameter on the simulation model. For example, when the number of parameters to be estimated is 15, the parameter space dimension is 15. This application uses Bayesian optimization to find the optimal parameter combination, and the number of simulation model calls is 10. 3 High-precision parameter estimation can be achieved at the order of magnitude. Compared with traditional optimization methods, this method significantly reduces the number of simulation model calls. BRIEF DESCRIPTION OF THE DRAWINGS

[0088] Figure 1 Schematic diagram of the optimization method of the present invention;

[0089] Figure 2 This is a key phase signal and a partial enlarged diagram of an embodiment of the present application;

[0090] Figure 3 is a transient vibration signal of an embodiment of the present application;

[0091] Figure 4 It is a comparison diagram of the simulated transient vibration curve and the measured vibration curve described in the embodiment of the present application; the dotted line is obtained by simulation, and the rotor dynamics simulation model is modeled using the optimal parameters obtained by this method; the solid line is the measured transient vibration signal. DETAILED DESCRIPTION

[0092] The technical solutions of the present application will be described clearly and completely below in conjunction with specific embodiments. However, those skilled in the art will understand that the embodiments described below are part of the embodiments of the present application, not all of them, and are only used to illustrate the present application and should not be considered to limit the scope of the present application. Based on the embodiments in the present application, all other embodiments obtained by those of ordinary skill in the art without creative work are within the scope of protection of the present application.

[0093] In one embodiment, referring to Figure 1 As shown, a Bayesian optimization estimation method for the dynamic parameters of a rotating machine is provided, comprising: synchronously collecting the key phase signal S when the rotating machine speed changes tacho And transient vibration signal x(t), and reconstruct the speed pulse time series and the time domain signal x clean (t); for vibration signal x clean (t) performing Fourier transform to obtain transient spectrum data, and extracting transient response data of a specified order therefrom; deleting the displacement offset caused by shaft bending and / or installation deviation in the transient response data; performing function fitting on the transient response data after eliminating the offset, and generating a transient response function g with the rotation speed or frequency as the independent variable and the amplitude as the dependent variable; i (ω); define the optimization parameter space and objective function; pass the objective function and optimization parameter space to the Bayesian optimization model to obtain the optimal parameter combination, and identify the rotating machinery dynamic parameters based on the optimal parameter combination. This method uses the transient vibration data of the rotating machinery and the simulated rotor model to construct the objective function. By applying Bayesian optimization to the parameter estimation of the rotating machinery rotor system, it solves the problem of dynamic parameter identification of complex rotor systems. It is implemented in the following steps:

[0094] 1. Collect the transient vibration signal x(t) and key phase signal S when the actual rotating machinery speed changes tacho , the corresponding time series t. Among them, the transient vibration signal x(t) and the key phase signal S tacho Synchronous acquisition: the signal changes obtained by the sensor are converted into electrical signals through the acquisition module, transmitted to the computer storage medium and saved.

[0095] Typical key phase signal waveform reference Figure 2 As shown in the figure, there are two main ways to obtain key phase signals:

[0096] First, the built-in hardware of the rotating machinery (such as an encoder) directly outputs a step-like pulse signal representing the rotor phase information. This hardware generates an electrical signal that strictly corresponds to the rotor's rotation period by synchronizing the mechanical structure with the rotor's motion.

[0097] Second, based on the periodic structural features (such as keyways and bosses) pre-set on the surface of the rotating machinery's shaft or turntable, a non-contact eddy current sensor or photoelectric sensor is used to collect the periodic step signals generated during the rotor's rotation. The generation of this signal depends on the relative position of the sensor and the rotor surface features.

[0098] Typical transient vibration signal waveform reference Figure 3 As shown in the figure, the transient vibration signal acquisition solution includes the following two technical paths:

[0099] Non-contact eddy-current sensors are placed on rotating machinery components, such as shafts and turntables, to measure the relative displacement of these components during operation. This technology uses the principle of electromagnetic induction to obtain real-time information on the relative displacement between the component being measured and the sensor probe, making it suitable for complex operating conditions such as high temperatures and high speeds.

[0100] For supporting structures such as bearing seats, contact accelerometers are used to measure vibration acceleration. These sensors convert mechanical vibrations into electrical signals through the piezoelectric effect, accurately capturing the dynamic response characteristics of the rotor system's supports. When the acquired signal is an acceleration signal, digital signal processing techniques are used to implement noise reduction filtering, and a quadratic integration algorithm is used to convert the acceleration signal into a displacement signal, enabling quantitative analysis of the vibration amplitude.

[0101] 2. Perform noise reduction on the key phase signal.

[0102] 2.1 Dynamic detection of the threshold value of the key phase pulse signal: Calculate the maximum value of the signal max(S tacho ) and the minimum value min(S tacho ), take the average value of peak and valley values ​​T threshold as the trigger threshold.

[0103]

[0104] 2.2 Capturing the leading edge of the key phase pulse signal: Using the double pointer sliding detection method to traverse the signal sequence, when S tacho (i)≤T threshold And S tacho (i+1)>T threshold When , record the index pair [i,i+1].

[0105] 2.3 Obtaining the pulse front time series of the key phase signal by interpolation Perform linear interpolation on each index pair that crosses the threshold to calculate the pulse start time

[0106]

[0107] Where Δt = 1 / fs , f s is the sampling frequency, and i represents the index in the sequence.

[0108] 2.4 According to the predefined speed threshold ω threshold , delete the low-speed signals at the beginning and end of the signal: by calculating the pulse front time series The first-order differential of is used to calculate the rotor system speed. Traverse the speed sequence from the beginning and end of the signal respectively, when it appears continuously When the index I is marked to be deleted del ; Delete by index and Pulse time points corresponding to the first and last low-speed zones and speed points Keep the valid data time series in the intermediate stable speed change area and speed sequence

[0109] 2.5 Double Smoothing Optimization of Time Series Speed ​​after cutting Apply Savitzky-Golay filter with window length N window,SG Set to 20, the polynomial order is 3, and the smooth speed series is obtained The i-th value in the sequence is:

[0110]

[0111] Among them, a j are the coefficients of the polynomial to be fitted, For sequence The jth value in .

[0112] 2.6 pairs Execute the filtering with the same parameters again to generate a secondary smoothed sequence

[0113] 2.7 Correction of time series based on smoothed speed: in, All are indexes, N3 is The sequence length.

[0114] 2.8 Output the speed pulse time series after noise reduction

[0115] 3. Perform noise reduction on vibration signals.

[0116] 3.1 Constructing the angular domain: Among them, N rev is the speed pulse time series Represents the total number of rotation cycles.

[0117] 3.2 Using cubic spline interpolation to generate angular domain time scale: The original data used for interpolation is The sequence to be interpolated is:

[0118]

[0119] in, f θ The set angular domain sampling rate is generally 100, and the interpolation is the same as The corresponding angular domain time series

[0120] 3.3 Based on angular domain time series The original transient vibration signal x(t) is resampled and the segmented Hermite interpolation is used during sampling. After interpolation, the angular domain time series is obtained. The corresponding angular domain vibration signal x(θ).

[0121] 3.4 Establish a sliding window synchronous averaging window model and define the window step size as: S step , usually set to 1.

[0122] 3.5 Perform sliding window synchronous averaging on the diagonal domain vibration signal x(θ): take N avg The signals of adjacent cycles are averaged in the angular domain to obtain a signal with a length of one cycle after noise reduction, and then the signal is averaged with a step size of S. step Slide the window and perform the next average to finally obtain the angular domain signal x after noise reduction clean (θ). Where N avg is the average number of cycles, usually 7. The cycle in this step is the angular domain cycle 2π.

[0123] 3.6. Since the signal length after noise reduction will be shortened, according to the angular domain signal x after noise reduction clean The length of (θ) generates the corresponding angular domain time series and time series in the time domain

[0124] 3.7. Reconstructing the time domain time series: Reconstructing the new angular domain time series Perform uniform sampling with the same sampling frequency as the transient vibration signal x(t) to obtain the time series t 5 .

[0125] 3.8. Restore the time domain signal: Use piecewise cubic spline interpolation to convert the angular domain signal x clean (θ) is mapped to an equally spaced time series t 5 , get the denoised time domain signal x clean(t); the original data used for interpolation is The sequence to be interpolated is t 5 .

[0126] 4. Extract the transient response of the specified order through Fourier transform.

[0127] 4.1 Transform the time domain signal x clean (t) Transform from time domain to angle domain signal x clean (θ), the sampling frequency in the angular domain is f s , the method is the same as steps 3.1 to 3.3.

[0128] 4.2 Determine the dynamic fundamental frequency: Based on the time series in the time domain Calculate the dynamic fundamental frequency period:

[0129]

[0130] Among them, N base is the set fundamental frequency denominator parameter, which is 1 in this embodiment.

[0131] 4.3 For each fundamental frequency period, the angular domain signal x clean (θ) performs a windowed fast Fourier transform (FFT) to obtain the complex matrix X:

[0132] [X=[X1,X2,…,X n ] T ,(n <f s / 2);

[0133] The i-th row in the matrix represents the rotation speed i / N base Times the frequency domain component of the frequency vibration.

[0134] 4.4 Extract the required rows from the complex matrix X and extract the magnitudes of all complex elements in the row |X i,j | and stored as |X| matrix:

[0135]

[0136] Each row of data in the matrix represents the rotation speed i / N that changes with the rotation phase base times the amplitude of the frequency vibration, N Xij The number of angular domain signal data points within each fundamental frequency period.

[0137] 4.5 Construct the angular velocity domain corresponding to the |X| matrix by Differential calculation of instantaneous speed ω(i):

[0138]

[0139] Then the instantaneous speed ω(i) is smoothed in the same way as steps 2.5 and 2.6 to obtain the smoothed instantaneous speed vector ω corresponding to the matrix |X| 2 .

[0140] 5. Remove the transient response bias in the |X| matrix to obtain |X 2 |:For each row X in |X| matrix i Remove the bias separately, the strategy is as follows:

[0141]

[0142] Among them, X i,j is the element in row i and column j of |X| matrix; For |X 2 | is the element in the i-th row and j-th column of the |X| matrix; β is the attenuation factor, usually 0.9; T is the adaptive threshold, usually the value of the first element in the i-th row of the |X| matrix.

[0143] 6. Transient response functionalization: Take |X 2 Each row of data in the matrix and the corresponding instantaneous speed vector ω 2 As the source data, the cubic spline interpolation method is used to construct the proxy function of the source experimental data In order to achieve smoothness during cubic spline interpolation, the following objective function is minimized and the transient response function g is stored. i (ω):

[0144]

[0145] Wherein, p is the smoothing parameter, which is 0.8 in this embodiment; w j is the weight function, which is 1 in this embodiment; is the second-order derivative of the surrogate function.

[0146] 7. Define the optimization parameter space.

[0147] 7.1. In this embodiment, the number of optimized parameters is 16, which are support stiffness k b1 ,k b2 ,k b3 ,k b4 , support damping c b1 ,c b2 ,c b3 ,c b4 , boundary stiffness k e1 ,k e2 ,k e3, turntable imbalances e1, e2, e3, e4, and structural damping parameter β. These parameters are all model parameters in the simulated rotor model. The simulated rotor model should be a proxy model / digital twin of the actual rotating machinery and can represent most of the dynamic characteristics of the actual rotating machinery.

[0148] 7.2. Define the value range of each optimization parameter: Support stiffness value range [10 6 ,10 9 ], the support damping value range is [10,3×10 4 ], the boundary stiffness range is [10 4 ,10 7 ], the range of turntable unbalance value is [10 -5 ,10 -2 ], the range of structural damping parameters is [10 -5 ,3×10 -4 ].

[0149] 7.3. Perform logarithmic transformation and normalization on each optimization parameter: Convert the scale of the optimization parameter space to a logarithmic scale and then perform normalization. The specific expression is as follows:

[0150]

[0151] Where h(x) represents the logarithmic and normalized function; x represents the optimization parameter.

[0152] 8. Define the optimization objective function.

[0153] 8.1. Objective function reads the optimization parameters and constant parameters of the simulation rotor model.

[0154] 8.2. Denormalize and delogarithmize the optimization parameters:

[0155]

[0156] 8.3. Call the simulated rotor model, update the optimization parameters in the rotor model, and re-establish the rotor dynamics differential equation:

[0157]

[0158] Where y represents the rotor system node displacements, M, K, C, G, and N represent the rotor system mass, stiffness, damping, gyroscope, and transient matrices, respectively, and F represents the external force. These system matrices and vectors are affected by changes in the optimization parameters, and each iteration results in an update of these matrices and vectors. The simulated rotor model used is a proxy model / digital twin of the actual rotating machinery and can represent most of its dynamic characteristics.

[0159] 8.4. Using the solver to solve the vibration response of the rotor model sim When solving, the changes in the rotation phase, rotation speed, and rotation acceleration of the simulation model are consistent with those during the experimental test, that is, the operating state of the simulation model is consistent with the operating state of the rotating machinery when the transient vibration signal is obtained.

[0160] 8.5. From the simulated vibration response y sim Extract the response data of some nodes The locations of these nodes are consistent with the locations of the measured sensors in the rotating machinery.

[0161] 8.6. Response Data Detrending was performed by subtracting the mean to remove the effects of gravity.

[0162] 8.7. Response Data Fourier transform is used to extract the transient response of the specified order and obtain |X sim |Matrix and ω sim Vector, the calculation method is the same as step 4.

[0163] 8.8. Define the weight function and use the linear combination of multi-peak Gaussian functions as the weight function:

[0164]

[0165] Among them, w(Ω) is the weight coefficient at the speed Ω, α i is the amplitude coefficient of the i-th Gaussian kernel, μ i is the central speed of the i-th Gaussian kernel, σ i The standard deviation of the i-th Gaussian kernel is used for weighted bandwidth control. The weight function is normalized using the weight coefficient.

[0166] 8.9. Define the objective function as:

[0167]

[0168] Ω={ω sim (1),ω sim (2), …};

[0169] in, is the scalar value of the objective function f(x); N component ,N signal are the set number of frequency domain components and the set number of signals, respectively. In this embodiment, N component =1,N signal =8, that is, only the 1-fold frequency component of the speed and 8 different transient vibration signals are considered. The 8 different transient vibration signals come from different measuring points and working conditions of the rotating machinery; Ω represents the matrix |X sim|All elements in the corresponding instantaneous speed vector; g i,j (Ω) represents the value of the transient response function of the i-th frequency domain component and the j-th signal at the speed Ω; |X sim | i,j (Ω) represents the i-th frequency domain component, the j-th signal corresponds to the instantaneous speed Ω |X sim |The value of the matrix; w i,j (Ω) represents the weight coefficient of the i-th frequency domain component and the j-th signal at the instantaneous speed Ω.

[0170] 8.10. Output the objective function value corresponding to the optimization parameters.

[0171] 9. Execute Bayesian optimization and output the calculated optimal parameters if the termination condition is met.

[0172] 9.1. Read the parameter space defined in step 7 Wherein x is the optimization parameter, X is the parameter space, and D represents the data set. In this embodiment, the number of elements in the data set D is 16.

[0173] 9.2. Set the objective function Obey the Gaussian process prior:

[0174] f(x)~GP(m(x),k(x,x′;λ,l));

[0175] Where m(x) is the mean function; λ is the output scale hyperparameter that controls the output amplitude, and l is the length scale hyperparameter that controls the smoothness of the function. The covariance kernel function k(x,x') can be calculated using the Matérn 5 / 2 kernel:

[0176]

[0177] In each optimization iteration, new hyperparameters are obtained based on the data set D, and then a new Gaussian process is obtained to approximate the unknown objective function f(x).

[0178] 9.3. During each iterative optimization, a new Gaussian process is used to obtain the posterior distribution, providing the predicted value and variance. The acquisition function is then used to select the next point to be evaluated x∈X with the maximum benefit. The acquisition function can use the expected improvement (EI):

[0179] α EI =E[max(f(x * )-f(x),0)]=∫ max (f(x * )-f(x),0)p(f(x)|x)df(x);

[0180] Among them, E[·] is the expectation, x *is the optimal set of parameters currently observed, and p(·) represents the probability distribution. The next evaluation point is selected by maximizing the acquisition function:

[0181]

[0182] 9.4. Select the evaluation point x next Substitute the objective function f(x) to calculate the scalar value L next , and the newly obtained data {x next ,L next Update to dataset D:

[0183] D′=D∪{x next ,L next}

[0184] After the dataset is updated, proceed to step 9.2.

[0185] 9.5. When the set termination condition is reached, stop the optimization iteration and output the optimal parameter combination x currently observed * The termination condition is usually the maximum number of evaluations, evaluation time limit, convergence threshold, etc.

[0186] This embodiment uses the maximum evaluation number of 5000 as the termination condition. The comparison between the simulated transient vibration curve corresponding to the optimal observation parameters and the measured vibration curve is shown in Figure 4 From the figure, we can see that through this method, the simulated rotor model can simulate most of the dynamic characteristics of the actual rotating machinery, proving that the dynamic parameters estimated by this method are accurate; the results in the figure are obtained after 5000 optimization evaluations. Considering that the parameter space dimension is 16, the results in the figure can prove that this method has high evaluation efficiency.

[0187] Although the embodiments of the present application are described above in conjunction with the accompanying drawings, the present application is not limited to the above-mentioned specific embodiments and application fields. The above-mentioned specific embodiments are merely illustrative and instructive, and not restrictive. A person of ordinary skill in the art, guided by this specification and without departing from the scope of protection of the claims of this application, may also devise various forms, all of which fall within the scope of protection of this application.

Claims

1. A Bayesian optimization estimation method for dynamic parameters of rotating machinery, characterized in that: include: S1: Synchronously collect the key phase signal S when the rotating machinery speed changes tacho And transient vibration signal x(t), and reconstruct the speed pulse time series and the time domain signal x clean (t); S2: vibration signal x clean (t) performing Fourier transform to obtain transient spectrum data, and extracting transient response data of a specified order therefrom; S3: deleting the displacement offset caused by shaft bending and / or installation deviation in the transient response data; S4: Perform function fitting on the transient response data after eliminating the bias, and generate a transient response function g with the speed or frequency as the independent variable and the amplitude as the dependent variable. i (ω); S5: Define the optimization parameter space and objective function: The optimization parameter space includes at least one of support stiffness, structural damping parameter, rotor eccentricity, bearing force parameter and oil film force parameter; Based on the transient response function g i (ω) Generate an objective function that receives the optimization parameters and calls the simulated rotor model to output a scalar value: Oh = {oh sim (1),oh sim (2),…}; in, is the scalar value of the objective function f(x); N component ,N signal are the set number of frequency domain components and the set number of signals respectively; Ω represents the matrix |X sim All elements in the corresponding instantaneous speed vector, matrix |X sim | is the complex matrix of the specified order transient response extracted from the simulated vibration response data by Fourier transform in the vibration response analysis of the rotor model. The elements in the matrix represent the amplitude of the corresponding order vibration component at different frequencies; g i,j (Ω) represents the value of the transient response function of the i-th frequency domain component and the j-th signal at the speed Ω; |X sim | i,j (Ω) represents the i-th frequency domain component, the j-th signal corresponds to the instantaneous speed Ω |X sim |The value of the matrix; w i,j (Ω) represents the weight coefficient of the i-th frequency domain component and the j-th signal at the instantaneous speed Ω; S6: passing the objective function and the optimization parameter space to the Bayesian optimization model to obtain the optimal parameter combination, and identifying the dynamic parameters of the rotating machinery based on the optimal parameter combination.

2. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Denoising and reconstructing speed pulse time series include: S111: Get key phase signal S tacho The average value of the peak and valley values ​​T threshold as a trigger threshold; S112: Get key phase signal S tacho Pulse front time series Calculate the pulse start time based on the trigger threshold Where Δt = 1 / f s , f s is the sampling frequency, i represents the index in the sequence; S113: According to the predefined speed threshold ω threshold Delete the low-speed signals at the beginning and end of the signal to generate a valid data time series and speed sequence S114: Speed ​​sequence Double smoothing optimization time series Generate a quadratically smoothed speed series S115: Based on smooth speed sequence Correct the time series to generate the speed pulse time series 3. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Denoise and reconstruct the time domain signal x clean (t) include: S121: Based on the speed pulse time series Constructing the corner area Use cubic spline interpolation to generate The corresponding angular domain time series S122: Based on angular domain time series Mapping the transient vibration signal x(t) to the angular domain vibration signal x(θ); S123: Use a sliding window to perform synchronous averaging on the angular domain vibration signal x(θ) to generate an angular domain signal x clean (θ); S124: Based on the angle domain signal x clean (θ) generates the corresponding angular domain time series and time series in the time domain S125: According to the frequency of the transient vibration signal x(t), the time series Uniform sampling generates time series t 5 ; S126: The angle domain signal x clean (θ) is mapped to the time series t 5 Generate time domain signal x clean (t).

4. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Step S2 includes: S21: Based on the speed pulse time series and the angular domain vibration signal x(θ) to generate the angular domain signal x clean (θ), and according to the angular domain signal x clean (θ) corresponds to the time domain time series Construct the fundamental frequency cycle: Among them, N base is the set fundamental frequency denominator parameter; S22: The time domain signal x clean (t) is mapped to the angular domain to generate the angular domain signal x clean (θ), and divided into blocks according to the fundamental frequency period; S23: For each block of inner angle signal x clean (θ) performs a windowed fast Fourier transform and outputs a complex matrix X: X=[X1,X2,…,X n ] T ,(n <f s / 2)? The i-th row in the matrix represents the rotation speed i / N base times the frequency domain component of the frequency vibration; S24: extract the complex elements of the row of the specified order from the complex matrix X, calculate the amplitude and generate an amplitude sequence; S25: Construct the angular velocity domain corresponding to the amplitude sequence, based on the time domain time series The instantaneous speed ω(i) is calculated by differential calculation and the instantaneous speed vector ω is generated by smoothing. 2 , forming the transient response data of amplitude-speed.

5. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Step S3 includes: extracting amplitude values ​​from the transient spectrum data and storing them as a |X| matrix; for each row X in the |X| matrix, i Remove the offsets individually: Among them, X i,j is the element in the i-th row and j-th column of the |X| matrix; β is the attenuation factor; and T is the adaptive threshold.

6. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Step S4 includes: S41: Obtain the row vector of the amplitude sequence and the corresponding instantaneous speed vector ω based on the transient response data 2 ; S42: With speed as the independent variable and amplitude as the dependent variable, a proxy function is constructed using cubic spline interpolation: S43: Minimize the following objective function and store the transient response function g i (ω): Where p represents the smoothing parameter; w j represents the weight function; represents the second-order derivative of the surrogate function; X i,j represents the element in the i-th row and j-th column of the amplitude sequence, i.e., the vibration amplitude corresponding to order i at the j-th instantaneous speed point; g i (ω) represents the predicted value of the surrogate function.

7. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Defining the optimization parameter space in step S5 includes: S51: defining optimization parameters and the value range of each optimization parameter; S52: Convert each parameter to the normalized space through the function h(x): Where h(x) represents the logarithmic and normalized function; x represents the optimization parameter.

8. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Defining the objective function in step S5 includes: S53: performing inverse normalization and inverse logarithmization on the optimized parameters to restore the physical parameters; S54: Input physical parameters into the rotor dynamics differential equation to update the rotor system mass matrix, stiffness matrix, damping matrix, gyro matrix, and transient matrix: Where y represents the node displacement of the rotor system; M, K, C, G, and N represent the mass matrix, stiffness matrix, damping matrix, gyro matrix, and transient matrix of the rotor system, respectively; and F represents the external force. S55: Determine the vibration response y under the same measured speed / acceleration conditions sim , extract the response data matching the measurement point S56: Response data Detrending is performed and Fourier transform is used to extract the transient response of the specified order to generate the simulation amplitude matrix |X sim |Matrix and speed vector ω sim ; S57: Generate an objective function using a multi-peak Gaussian function weighting.

9. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 8, characterized in that: The multi-peak Gaussian weight function is defined as: Among them, w(Ω) is the weight coefficient at the speed Ω, α i is the amplitude coefficient of the i-th Gaussian kernel, μ i is the central speed of the i-th Gaussian kernel, σ i is the standard deviation of the i-th Gaussian kernel used for weighted bandwidth control.

10. The Bayesian optimization estimation method for rotating machinery dynamic parameters according to claim 1, characterized in that: Step S6 includes: S61: Read the optimization parameter space and calculate the objective function Establish a Gaussian process prior model: f(x)~GP(m(x),k(x,x′;λ,l)) Where m(x) is the mean function; λ is the output scale hyperparameter that controls the output amplitude; l is the length scale hyperparameter that controls the smoothness of the function; The covariance kernel function k(x,x') is calculated using the Matérn 5 / 2 kernel function: Where D represents the dataset; S62: Calculate the Gaussian process posterior distribution based on the current observation data set D, and select the next evaluation point x by the expected improvement acquisition function next : α EI =E[max(f(x * )-f(x),0)]=∫max(f(x * )-f(x),0)p(f(x)|x)df(x) Among them, E[·] is the expectation, x * is the optimal set of parameters currently observed, and p(·) represents the probability distribution; S63: x next Substitute the objective function f(x) to calculate the scalar value L next , and the new data {x next ,L next }Update to dataset D to form a new dataset D`; S64: When the maximum number of evaluations or the convergence threshold is met, the current optimal parameter x is output. * .

Images