A Multivariate Feature Fusion Method for Robot Milling Edge Flutter Based on PSO-SVM-RFE

Abstract

The present invention discloses a multi - feature fusion method for robot milling edge chatter based on PSO - SVM - RFE. First, the acceleration signal is pre - processed to highlight the initial chatter frequency by eliminating the tooth frequency and its multiples. Secondly, the pre - processed acceleration signal is subjected to time - domain, frequency - domain and time - frequency domain transformations to obtain a feature vector set composed of waveform factor, center frequency and energy entropy. Then, based on the support vector machine recursive feature elimination method (SVM - RFE), the sensitivity of different feature quantities to milling edge chatter is studied to construct a multi - feature classification model. Next, the particle swarm optimization algorithm (PSO) is used to fuse the top 3 sensitive features before sorting into a new chatter feature, and the milling edge process is divided into three stages: stable state, initial chatter and severe chatter by the chatter threshold. Finally, the above algorithm program is written to calculate the fusion features, and the recognition results are verified through the experimental data of robot milling edge. The present invention can effectively solve the problem of feature quantity selection of robot milling edge chatter signals and improve the chatter recognition accuracy.

Description

Technical Field

[0001] The present invention belongs to the technical field of robot edge milling, and specifically relates to a multi - feature fusion method for robot edge milling chatter based on PSO - SVM - RFE. Background Technique

[0002] Currently, carbon fiber reinforced polymer (CFRP) has been favored in the fields of aerospace, military, and automotive due to its excellent properties such as designability, light weight, high strength, fatigue resistance, corrosion resistance, and easy repair. However, due to its poor forming accuracy, a margin is usually left as a sacrificial layer at the edge of the component during the forming process, and it is removed by milling after curing, which is called edge milling, to meet the assembly and connection requirements with other components. For the edge milling of CFRP parts, the existing traditional machine tools have limitations in the processing stroke, the cost of developing special machine tools is high, and the manual edge cutting has low efficiency and poor accuracy. Therefore, it has become a trend to use industrial robots with high flexibility and low cost for large - scale composite material edge milling. However, the disadvantage of its weak structural rigidity leads to the easy occurrence of chatter during the edge milling process. Strong vibrations not only damage the quality of the machined surface, affect the edge milling accuracy of the parts, but also exacerbate tool wear and damage the robot edge milling system. In recent years, using sensors to collect signals related to the robot processing state, extracting characteristic quantities reflecting chatter information, and realizing real - time monitoring of robot processing chatter based on artificial intelligence algorithms have gradually become new concerns of domestic and foreign scholars. However, too many chatter characteristic quantities will produce information redundancy during the realization of chatter monitoring, which will affect the recognition accuracy and duration at the same time. Using a feature selection algorithm to screen out the features sensitive to chatter and fusing them into one feature can improve the chatter recognition accuracy rate and also contribute to the determination of the threshold. This fused feature has the advantages of multi - features on the one hand. Compared with single features, it has strong anti - interference ability and high accuracy. On the other hand, it avoids the situation that when using multi - features to monitor chatter, it can only reflect whether there is chatter or not, but cannot reflect the severity of chatter. By monitoring the comparison between the value of the fused feature and the threshold, the severity of chatter can be reflected, providing a reference for future chatter suppression. Summary of the Invention

[0003] The purpose of the present invention is to provide a multi - feature fusion method for robot edge milling chatter based on PSO - SVM - RFE, which combines machine learning algorithms with chatter monitoring to optimize the selection of edge milling chatter characteristic quantities and improve the recognition accuracy rate of robot milling chatter.

[0004] The technical solution to achieve the above purpose is as follows:

[0005] A multi - feature fusion method for robot edge milling chatter based on PSO - SVM - RFE includes the following steps:

[0006] Step 1: Preprocess the vibration signals collected by the acceleration sensor: Use a Butterworth filter to remove high-frequency noise signals; perform spectral analysis on the denoised signals, and use the Frequency Elimination Algorithm (FEA) to remove the tooth passing frequency and its multiples.

[0007] Step 2: Extract the chatter characteristics of the acceleration signal: Perform time-domain, frequency-domain, and time-frequency-domain transformations on the preprocessed acceleration signal respectively to extract the corresponding time-domain, frequency-domain, and time-frequency-domain characteristic quantities.

[0008] Step 3: Rank the weights of the feature vectors based on SVM-RFE: Based on the large number of characteristic quantities obtained in Step 2, combine them into a feature vector set, and use the Support Vector Machine Recursive Feature Elimination (SVM-RFE) method. Take the feature quantity vector set as the input to obtain the sorted feature vector set.

[0009] Step 4: Multivariate feature fusion: Based on the Particle Swarm Optimization Algorithm (PSO), perform multivariate feature fusion on the top three sensitive features in the feature vector set obtained in Step 3, with the maximum cost function as the optimization goal. Substitute the fused features into the SVM for verification of the chatter monitoring accuracy.

[0010] Compared with the prior art, the advantages of the present invention include the following four points:

[0011] (1) The present invention proposes to combine machine learning algorithms with chatter monitoring, optimizes the selection problem of chatter characteristic quantities, and has a good effect on chatter recognition.

[0012] (2) This method realizes the recognition of the severity of chatter and successfully solves the problem that the traditional method can only recognize whether there is chatter or not.

[0013] (3) A method for removing the tooth passing frequency and its multiples is proposed to solve the problem that it is difficult to monitor the initial stage of chatter, and the accuracy of monitoring the initial stage of chatter is improved.

[0014] (4) A method of multivariate fusion features is proposed to solve the problem that single feature recognition of chatter is easily interfered. It not only improves the recognition accuracy but also enhances the anti-interference ability.

[0015] In addition to the objectives, features, and advantages described above, the present invention has other objectives, features, and advantages. The following will refer to the accompanying drawings to further elaborate on the present invention in detail. Description of the Drawings

[0016] Figure 1 It is a flowchart of the stability analysis method of the present invention.

[0017] Figure 2 It is a diagram of the FEA result of the simulation signal.

[0018] Figure 3 It is the time-domain feature diagram of the typical chatter signal of robot edge milling.

[0019] Figure 4 It is the frequency-domain feature diagram of the typical chatter signal of robot edge milling.

[0020] Figure 5 It is the diagram of the correlation coefficient and energy ratio variance of each frequency band after FEA.

[0021] Figure 6 It is the time-frequency domain feature diagram of the typical chatter signal of robot edge milling.

[0022] Figure 7 It is the schematic diagram of the optimal decision surface of SVM.

[0023] Figure 8 It is the accuracy diagram of the feature subset model.

[0024] Figure 9 It is the diagram of the typical chatter signal of robot milling.

[0025] Figure 10 It is the spectrogram of the typical chatter signal of robot milling.

[0026] Figure 11 It is the time-domain and frequency-domain diagram of the typical chatter signal after FEA.

[0027] Figure 12 It is the change curve of the fusion feature a17. Specific implementation manner

[0028] The present invention will be further described below with reference to the accompanying drawings of the specification.

[0029] Figure 1 It is the flowchart of the modeling process of the prediction method of the present invention, specifically including the following steps:

[0030] Step 1. Preprocess the vibration signal collected by the acceleration sensor: Use a Butterworth filter to remove high-frequency noise signals; perform spectral analysis on the denoised signal, and use the frequency removal algorithm (FEA) to remove the tooth passing frequency and its multiples;

[0031] Step 1.1. Use a Butterworth filter to remove high-frequency noise signals. The Butterworth low-pass filter can be expressed by the following formula of the square of the amplitude with respect to the frequency:

[0032]

[0033] where n = the order of the filter, ω c is the cut-off frequency, ω pis the edge frequency of the passband. Since the flutter frequency is below 4000 Hz, a filter is used to filter out the high-frequency noise signal greater than 4000 Hz;

[0034] Step 1.2: Perform FFT on the acceleration signal {s(k), k = 1, 2,..., M} after removing the high-frequency noise signal to obtain the spectrum S(i), as shown in Equation (2). The length of the spectrum is M = length(X):

[0035]

[0036] where j is the imaginary unit;

[0037] Step 1.3: Set the tooth passing frequency f z , and its multiple frequency is c·f z , with the initial value c = 1. As shown in Equation (3):

[0038] f z = Ω / 60 (3)

[0039] where Ω is the spindle speed;

[0040] Step 1.4: Considering that the tooth passing frequency f z fluctuates due to other factors such as machining load, introduce a floating constant a, with the default value a = 0.5. Set the spectrum amplitude in the interval corresponding to the frequency (f z -a, f z +a) to zero, as shown in Equation (4);

[0041]

[0042] Step 1.5: If i is greater than c·f z +a and less than N / 2, where N is the number of signal points after Fourier transform, then let c = c + 1, and repeat Steps 1.3 to 1.5; otherwise, proceed to the next step;

[0043] Step 1.6: Since the spectrum is symmetric about the center after FFT, let to obtain the complete spectrum after removing the tooth passing frequency and its multiple frequencies, as Figure 2 shown;

[0044] Step 1.7: Perform inverse Fourier transform (IFFT) on the spectrum after removing the tooth passing frequency and its multiple frequencies to obtain the frequency-removed signal s(k), as shown in Equation (5):

[0045]

[0046] Step 2. Flutter feature extraction of acceleration signal: Perform time-domain, frequency-domain, and time-frequency domain transformations on the preprocessed acceleration signal respectively to extract the corresponding time-domain, frequency-domain, and time-frequency domain feature quantities;

[0047] Step 2.1. Extract the time-domain features of the acceleration signal after frequency removal: It is found in the variable cutting depth experiment that when the cutting depth gradually increases, it can be clearly seen in the time domain that the signal amplitude continuously increases. However, it is difficult to accurately judge the robot milling edge state only based on the signal amplitude feature. Usually, dimensional or dimensionless parameters are extracted from the time domain as flutter monitoring features. Let the collected signal time series be {s k (t), k = 1, 2, …, M}, and the extracted time-domain feature set be T d , and the features are as follows:

[0048] (1) Maximum value a 1 : Describes the highest peak value of the vibration signal during the milling process, as shown in Equation (6):

[0049] a 1 = max(|s k |) (6)

[0050] (2) Average value a 2 : Describes the average change trend of the vibration signal during the milling process, as shown in Equation (7):

[0051]

[0052] (3) Peak-to-peak value a 3 : Describes the difference between the highest peak and the lowest peak of the vibration signal during the milling process, as shown in Equation (8):

[0053] a 3 = a 1 - min(s k ) (8)

[0054] (4) Variance a 4 : Describes the severity of the amplitude fluctuation change of the milling vibration signal, as shown in Equation (9):

[0055]

[0056] (5) Standard deviation a 5 : Describes the degree of dispersion between signal points of the milling vibration signal, as shown in Equation (10):

[0057]

[0058] (6) Kurtosis a 6: Describes the kurtosis of the probability distribution of vibration signal points. It is a dimensionless parameter, independent of the amplitude of the signal, but sensitive to the impact characteristics of the signal. a 6 = 3 indicates that the peak value of the signal distribution is close to that of the normal distribution at this time; when a 6 > 3, it means that the peak value of the signal distribution is greater than that of the normal distribution at this time, that is, the degree of dispersion of the signal points decreases, also known as positive kurtosis; when a 6 < 3, it is negative kurtosis. At this time, the degree of dispersion of the signal distribution increases, and the peak value of the signal distribution is lower than the peak value of the normal distribution curve, as shown in Equation (11):

[0059]

[0060] (7) Skewness a 7 : Describes the degree of asymmetry of the distribution of milling vibration signal points. Corresponding to kurtosis, kurtosis describes the situation of the peak value of the signal point distribution and the normal distribution, while skewness describes the degree of symmetry of the tails on both sides of the peak value of the signal distribution. Therefore, it can be known that in the case of normal distribution, a 7 = 0, that is, the two sides of the peak value are completely symmetric. When a 7 < 0, it is called left skewness, that is, at this time the signal points are mainly located on the right side of the mean value, showing that the left tail of the signal point distribution curve is longer; when a 7 > 0, it is called right skewness, and at this time the tail on the right side of the distribution curve is longer. Skewness is also a dimensionless characteristic, as shown in Equation (12):

[0061]

[0062] (8) Root mean square a 8 : Describes the average power and average energy of the vibration signal. Different from the mean value, the root mean square value first squares and sums the data, then averages, and finally takes the square root. It can not only describe the situation of the mean value to a certain extent, but also reflect the degree of data dispersion, as shown in Equation (13):

[0063]

[0064] (9) Form factor a 9 : Describes the waveform change trend of the milling vibration signal, as shown in Equation (14):

[0065]

[0066] (10) Crest factor a 10 : Describes the vibration moment situation of the milling vibration signal, as shown in Equation (15):

[0067]

[0068] (11) Impulse factor a 11: Describe the sensitivity of the milling vibration signal to pulses, as shown in Equation (16):

[0069]

[0070] (12) Margin factor a 12 : Describe the impact on the milling vibration signal, as shown in Equation (17):

[0071]

[0072] The above are the time-domain features extracted from the signal time domain, as Figure 3 shown, which include dimensional and dimensionless features. Therefore, the initial time-domain feature set T d is:[[]]

[0073] T d ={a 1 , a 2 ,..., a 12} (18)

[0074] Step 2.2: Extract the frequency-domain features of the acceleration signal, as Figure 4 shown: Different from the time-domain analysis that reflects the change trend of the signal amplitude over time, the frequency-domain analysis reflects the frequency information contained in the signal from another dimension. When chatter occurs during robot milling, not only will the signal time domain change significantly, but chatter frequencies close to the natural frequency of the robot system will also be generated in the frequency spectrum. The amplitudes of these frequencies directly indicate the instability degree of the system at this time. Therefore, extracting the spectral features in the signal is an important means for robot milling chatter detection. FFT is the main method for frequency-domain analysis, and its transformation formula is shown in Equation (2). After the signal is FFT-transformed, the frequency vector is {f i , i = 1, 2,..., M}, where f i represents the i-th frequency value of the frequency vector, and the amplitude of the frequency power spectrum corresponding to this frequency is b i , that is, the abscissa and ordinate in the spectrogram. Let the frequency-domain feature set composed of frequency feature quantities be T f , which is specifically as follows:

[0075] (1) Centroid frequency a 13 : Describe the frequency centroid in the frequency spectrum of the robot milling signal, that is, reflect the frequency distribution in the spectrum at a certain moment of the signal. In the case of stable robot milling, the signal spectrum mainly contains the tooth passing frequency and its multiples. When chatter occurs in the system, chatter frequencies will be generated in the spectrum. In the early stage of chatter, the frequency amplitudes are very small and are often submerged in the periodic signal. Therefore, the FEA algorithm can effectively remove the relevant harmonic frequencies. When the system generates chatter frequencies, the prominent chatter frequencies will inevitably cause the frequency centroid to shift. The centroid frequency formula is shown in Equation (19):

[0076]

[0077] (2) Mean square frequency a 14 : It describes the change of the main frequency band of the power spectrum in the robot milling signal spectrum. It has the same function as the centroid frequency, as shown in Equation (20).

[0078]

[0079] (3) Frequency variance a 15 : It describes the discreteness of the energy distribution of the power spectrum of the robot milling signal. It can be seen from the formula that the frequency variance is closely related to the centroid frequency, as shown in Equation (21):

[0080]

[0081] Frequency domain analysis transforms the continuous time-domain signal into the discrete frequency domain through FFT, which can more accurately analyze the content contained in the signal. The frequency domain feature set T f is as follows:

[0082] T f = {a 13 , a 14 , a 15} (22)

[0083] Step 2.3, perform wavelet packet decomposition on the acceleration signal after frequency removal: Use wavelet packet decomposition (WPD) to decompose the signal into multiple sub-signals, and use energy entropy as the time-frequency domain feature. When using WPD to perform multi-layer decomposition on the signal, the original signal will first be decomposed by the filter into an approximate signal and a detail signal, that is, the signal passes through the low-pass filter to obtain the low-frequency signal and through the high-pass filter to obtain the high-frequency signal. Then, at the next layer, the approximate signal and the detail signal need to be further decomposed, and the signal is divided into two parts: low frequency and high frequency again. WPD is shown in Equation (23), and the reconstruction of each node signal is shown in Equation (24):

[0084]

[0085]

[0086] In the formula, a is the high-pass filter coefficient, b is the low-pass filter coefficient, d is the wavelet packet decomposition coefficient, n is the decomposition layer number, and j is the frequency band node number. l = k + 1, k = M / 2 n , indicating the number of wavelet packet coefficients of each frequency band in the nth layer, and M is the number of sampling points of the time series S.

[0087] Step 2.4, Flutter-sensitive frequency band selection: To improve the recognition ability of the energy entropy of time-frequency domain feature quantities for the initial flutter frequency, multiple sub-signals after WPD are screened to remove the sub-signal frequency bands that are weakly sensitive to the flutter state. The correlation coefficient ρ and the variance σ of the fluctuation of the frequency band energy ratio are used i (t) 2 Select the sub-signals with strong sensitivity to robot milling flutter after WPD:

[0088] (1) Correlation coefficient

[0089] After the collected original acceleration signal is subjected to FEA and denoising, the main component contained in the de-frequenced signal is the flutter signal. The correlation coefficient ρ is used to reflect the correlation degree between each sub-signal s i and the original signal S, and its value range is [-1, 1], as shown in Equation (25).

[0090]

[0091]

[0092] where s i is the i-th sub-signal, S is the original signal, and cov(s i , S) is the covariance of s i and S. and are the means of s i and S respectively, D(s i ) and D(S) are the variances of s i and S respectively, and n is the number of wavelet packet decomposition layers.

[0093] Generally, it is considered that when |ρ(s i , S)| ≤ 0.1, it indicates that the sub-signal is slightly correlated with the original signal; 0.1 < |ρ(s i , S)| ≤ 0.5 is a real correlation; 0.5 < |ρ(s i , S)| ≤ 1 is a significant correlation. For the robot milling acceleration signal, the larger the absolute value of the correlation coefficient between the reconstructed signal of the sub-signal wavelet packet coefficient and the original signal, the greater the correlation degree and the more useful signal components it contains; on the contrary, it indicates that the signal contained in the frequency band is mainly interference noise.

[0094] (2) Variance of the fluctuation of the frequency band energy ratio

[0095] Under the stable milling state of the robot, the fluctuation of the energy ratio of each sub-signal band within each sampling period T follows a normal distribution, and the energy ratio of the sub-signal fluctuates up and down around the mean value. However, when chatter occurs during robot milling, the energy of some sub-signals will change sharply and finally gather at the natural frequency of the robot. The energy ratio of the sub-signal containing the chatter frequency increases, and the variance of the reconstructed signal of the wavelet packet coefficient of the sub-signal increases significantly in terms of the numerical value. Therefore, the variance of the energy ratio fluctuation σ i (t) 2 is used to represent the sensitivity of the reconstructed signal of the sub-signal to the abnormal chatter state, as shown in Equation (27):

[0096]

[0097] In the formula, μ i is the mean value of the energy ratio of the i-th sub-signal, and P i (t) is the energy ratio value of the i-th sub-signal at time t. The larger the variance σ i (t) 2 , the better the sensitivity of this sub-signal to robot milling chatter; on the contrary, it indicates a sub-signal with weak sensitivity to chatter. The calculated variances of the energy ratios of each sub-signal are rearranged in descending order as {σ 2} = {σ 1 2 , σ 2 2 , …, σ j 2 , j = 1, 2, …, 2 n}.

[0098] (3) Definition of chatter-sensitive frequency band

[0099] According to the correlation degree between each sub-signal and the original signal and the fluctuation of the energy ratio of the sub-signal with the chatter state, as Figure 5 shown, the sub-signals of the weak-sensitive frequency band are removed to reduce the redundant interference of the system calculation and improve the sensitivity of the time-frequency domain characteristic quantity to the monitoring of the initial frequency of robot milling chatter. Therefore, it is set that the chatter-sensitive sub-signal needs to meet the following two conditions simultaneously:

[0100] 1) The correlation coefficient is greater than [ρ]

[0101] ρ > [ρ] (28)

[0102] 2) The variance of the frequency band energy ratio fluctuation is among the top N i numerical values

[0103] {σ 2} = {σ 1 2 , σ 2 2 ,..., σ j2 , j = 1, 2, ..., N i}], N i ≤ 2 n (29)

[0104] N i = round(λ·N) (30)

[0105] where round is the rounding function in MATLAB and λ is the sensitivity coefficient.

[0106] Step 2.5, Extract the energy entropy of time-frequency domain feature quantities: Calculate the energy entropy for the selected sensitive frequency bands. The calculation formula of the energy entropy is shown in Equations (31) and (32):

[0107]

[0108]

[0109] In the formula, E i is the energy of the i-th frequency band, P i is the proportion of the energy of the i-th sub-signal after normalization in the total energy, and there is Set the correlation coefficient to [ρ] = 0.1, the sensitivity coefficient λ = 0.7, and N i is 11 according to Formulas (28) and (29). Use the "db3" wavelet basis to perform four-layer wavelet packet decomposition on the signal to obtain 16 sub-signals. According to the sensitivity coefficient and the correlation coefficient, finally eliminate the sub-signals s2, s9, s10, s11, and s12 and do not consider them when calculating the energy entropy of time-frequency domain features, as Figure 6 shown.

[0110] Step 3, Rank the weights of feature vectors based on SVM-RFE: Based on the large number of feature quantities obtained in Step 2, combine them into a feature vector set, and use the support vector machine recursive feature elimination method (SVM-RFE) to take the feature quantity vector set as the input to obtain the ranked feature vector set;

[0111] Step 3.1, Establish a support vector machine (SVM) classification model: As a binary classification machine algorithm, the main principle of SVM is to establish a classification hyperplane as the decision surface according to the distribution of data points, so that the isolation edge distance between the data points on both sides of the decision surface is maximized, as Figure 7 shown. Assume that the data set {(x i , y i ), i = 1, 2, …, n} is a set of training data, where i is the number of samples, x i is an n-dimensional feature vector, y iThe value range of the class label is {-1, 1}. The distance from the isolation edge L1 or L2 of the data point to the decision surface is defined as 1 / ||w||. Since the two boundaries are symmetric about the decision surface, the margin between L1 and L2 is 2 / ||w||. Therefore, the problem of finding the optimal separation plane can be transformed into solving for the minimum ||w||. The SVM decision function is shown in formula (33):

[0112] f(x i ) = ω·x i +b (33)

[0113] where ω = {ω 1 ; ω 2 ; …; ω d} is the weight vector, which determines the direction of the classification hyperplane. The d in the vector represents the number of eigenvalue. b is the displacement term, and the distance between the classification hyperplane and the origin is b / ||w||.

[0114] The position of the classification hyperplane depends on the weight vector ω and the displacement b. Considering that there are not always linearly separable data samples, the slack variable ξ i is introduced. By solving the following optimization constraint equation, the optimal classification hyperplane can be obtained.

[0115]

[0116] where ξ i = max(0, 1 - y i (ω·x i +b)), which allows the samples to not satisfy the constraints to a certain extent. Each data sample has a corresponding slack variable, indicating the degree to which the sample does not satisfy the constraints. C is the penalty function, and its value is greater than zero, representing the trade-off between the training error and the margin. The larger the C value, the greater the penalty for classification, the higher the fitting degree of the training samples, but the generalization ability of the model will be reduced. is the mapping function.

[0117] Using the Lagrange function to solve equation (34), the optimization constraint equation can be written as the following formula:

[0118]

[0119] where α i is called the Lagrange multiplier, which is a constant determined during the optimization process. To avoid the high dimensionality of the feature space and reduce the complexity of the classification model, the kernel function is introduced, which satisfies the Mercer theory and enables the training model to have a minimized algorithm.

[0120] Therefore, the SVM non-linear decision function can be described by formula (36):

[0121]

[0122] Among them, is the optimal vector obtained from Equation (35), In actual model training, the kernel function often selects the radial basis function, as shown in Equation (37):

[0123]

[0124] In the formula, σ is the width parameter of the Gaussian kernel function, which controls the radial action range of the kernel function. The generalization ability of the training model first increases and then decreases with the increase of σ. The iterative method is used to find the appropriate parameter σ to make the training model achieve the optimal generalization ability.

[0125] Step 3.2: Establish a support vector machine recursive feature elimination (SVM-RFE) model for feature vector sorting: The SVM-RFE algorithm uses the weight size of the SVM classifier as the feature sorting criterion. Based on the trained SVM classifier, the feature with the smallest weight is removed, and the remaining features are used as the input of the training model for a new round of iteration until all features are eliminated. The feature sorting result is obtained according to the elimination order. The cost function is used as the weight of each feature of the SVM. For a linear SVM, the cost function for eliminating the i-th feature is as shown in Equation (38):

[0126] DJ(i) = (ω i ) 2 (38)

[0127]

[0128] For a non-linear SVM, the cost function for eliminating the i-th feature is as shown in Equation (40):

[0129]

[0130] In the formula, DJ(i) is the cost function of the i-th removed feature. H is the matrix of elements y i y j K(x i , x j ). H(-i) is the matrix after eliminating the i-th feature. K is the kernel function, which represents the similarity between x i , x j .

[0131] Step 3.3: Perform flutter-sensitive feature selection for the robot milling edge based on the SVM-RFE model: Use the feature vector set composed of the 12 time-domain features, 3 frequency-domain features, and 1 time-frequency domain feature extracted in Step 2 as the input of the SVM-RFE. Before performing the SVM-RFE algorithm, first normalize the data to improve the calculation accuracy and time of the model. Use min-max normalization, as shown in formula (41):

[0132]

[0133] Then use the cross-validation algorithm (K-CV) to select the best penalty function and kernel function parameters. The penalty function C = 3.0314, and the Gaussian kernel function g = 4. After the SVM-RFE algorithm, the feature ranking vector r = [a 16 , a 9 , a 13 , a 14 , a 15 , a 6 , a 4 , a 8 , a 5 , a 3 , a 10 , a 7 , a 12 , a 2 , a 11 , a 1 is obtained.

[0134] Step 4: Multivariate feature fusion: Based on the particle swarm optimization algorithm (PSO), perform multivariate feature fusion on the top three sensitive features in the feature vector set obtained in Step 3, with the maximum cost function as the optimization goal. Substitute the fused features into the SVM for flutter monitoring accuracy verification.

[0135] Step 4.1: Select the features for fusion: Construct 16 feature subsets from the feature ranking vector obtained in Step 3 in sequence. Each subset adds one more feature compared to the previous subset. Substitute the 16 feature subsets into the SVM for classification training respectively. Finally, it is found that when the feature vector is energy entropy a 16 , centroid frequency a 13 , and waveform factor a 9 , the model accuracy is the highest, reaching 98.02%, and the training time is only 1 / 3 of that with full features, as Figure 8 shown.

[0136] Step 4.2: Perform multivariate feature fusion based on the particle swarm algorithm (PSO): Before iteration, the PSO algorithm initializes a group of random particles for the solution of the function, and then continuously iterates in the D-dimensional space to optimize the optimal position, that is, search for the global extreme value P g. During the iteration process, the particle continuously updates its velocity V and position X, and its "flight" change formula is as follows:

[0137]

[0138]

[0139] In the formula, ε is the inertia weight; i = 1, 2,..., n; k is the current iteration number; V in is the particle velocity; c 1 , c 2 is the acceleration factor, and its value range is [0, 2]. c 1 determines the particle's own best step size, and c 2 determines the particle's global best step size; r 1 , r 2 is a random number, and its value range is [0, 1].

[0140] Using the PSO algorithm, the sensitive feature set [a 16 , a 9 , a 13 selected in step 4.1 is fused into a fusion feature a17 that is strongly related to the flutter state. The fusion objective formula is as shown in formula (44):

[0141]

[0142] In the formula, a 17 is the new feature obtained by fusing features a 16 , a 9 and a 13 with certain coefficients. The sum of coefficients l 1 , l 2 , l 3 is 1. The optimization objective is to maximize the cost function of the fused feature, that is, to make the feature ranking of the fused feature be in the forefront after SVM - RFE.

[0143] Before optimizing the PSO - SVM - RFE algorithm, initialize c 1 = c 2 = 1.49445, k max = 100, n = 20. After optimization, the fusion coefficients are l 1 = 0.8393, l 2 = 0.1566, l 3 = 0.0041. The fused feature enables the accuracy of the SVM classifier to reach 99.25%, greatly improving the accuracy of the classifier.

[0144] The present invention sorts the chatter feature quantities of robot edge milling, and uses a fusion algorithm to fuse the feature quantities in the front row of the sorting into a new chatter monitoring feature quantity. The edge milling process is divided into three stages: a stable state, the initial stage of chatter, and severe chatter by the chatter threshold, realizing the accurate identification of the chatter state. It has stronger anti-interference ability compared with single feature quantities, and improves the accuracy of chatter prediction in robot edge milling processing, providing technical support and theoretical guidance for actual production.

[0145] Example 1

[0146] In this example, carbon fiber reinforced composite material (CFRP) is taken as the research object of edge milling. The model of the industrial robot is KUKAKR210 - R2700 EXTRA. The structural schematic diagram of the robot longitudinal-torsional ultrasonic edge milling system built is as Figure 7 shown. Its specific composition includes a robot body, a small longitudinal-torsional ultrasonic edge milling system, an ultrasonic generating device, and an end effector. The tool used in this research is a polycrystalline diamond (PCD) milling cutter, and the specific parameters are shown in Table 1. The experimental parameters of robot edge milling are shown in Table 2.

[0147] Table 1 Related parameters of PCD milling cutter

[0148]

[0149] Table 2 Robot chatter signal experiment

[0150]

[0151] According to the experimental parameters in Table 2, the experiment number (1) is repeated 4 times, and the experiment number (2) is repeated 3 times. The first three experimental data of each processing parameter are used as the SVM - RFE training data, and the fourth experiment of the experiment number (1) is used as the classifier model verification data. The specific results are as Figure 9 、 Figure 10 and Figure 11 shown. Among them Figure 9 represents the time domain diagram of the typical robot edge milling signal, which is divided into two states: stable and chatter. Figure 10 represents the frequency domain diagram of the typical robot edge milling signal, and the chatter frequency appears. Figure 11 is the time domain and frequency domain diagrams of the typical chatter signal after FEA. Compared with the original signal, the acceleration amplitude decreases and the chatter frequency becomes prominent. The time domain, frequency domain, and time - frequency domain features are respectively extracted from the de - frequencyed signal, as Figure 3 、 Figure 4 and Figure 6 shown. The fusion algorithm is used for the feature quantities after SVM - PFE to obtain the fusion feature a17, as Figure 12As shown, the threshold is set to 0.2, and the classification accuracy of the model can reach 99.25%, greatly improving the accuracy. Therefore, combining the PSO-SVM-RFE algorithm with the multi-features of robot milling edge chatter is a very effective way to solve the problem of robot milling edge chatter recognition.

[0152] The above are only the preferred embodiments of the present invention and are not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and changes. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principle of the present invention shall be included within the protection scope of the present invention.

Claims

1. A multi - feature fusion method for robot milling edge chatter based on PSO - SVM - RFE, characterized in that, it includes the following steps: Step 1: Pre - process the vibration signal collected by the acceleration sensor: Use a Butterworth filter to remove high - frequency noise signals; perform spectral analysis on the denoised signal, and use the frequency elimination algorithm (FEA) to remove the tooth passing frequency and its multiples; Step 2: Extract chatter features of the acceleration signal: Perform time - domain, frequency - domain, and time - frequency - domain transformations on the pre - processed acceleration signal respectively to extract the corresponding time - domain, frequency - domain, and time - frequency - domain feature quantities; Step 3: Rank the weights of feature vectors based on SVM - RFE: Based on the time - domain, frequency - domain, and time - frequency - domain feature quantities obtained in Step 2, combine them into a feature vector set, and use the support vector machine recursive feature elimination method (SVM - RFE). Take the feature quantity vector set as the input to obtain the ranked feature vector set; Step 4: Multi - feature fusion: Based on the particle swarm optimization algorithm (PSO), perform multi - feature fusion on the top three sensitive features in the feature vector set obtained in Step 3, with the maximum cost function as the optimization goal. Substitute the fused features into the SVM for chatter monitoring accuracy verification.

2. The multi - feature fusion method for robot milling edge chatter based on PSO - SVM - RFE according to claim 1, characterized in that, the specific steps of Step 1 are as follows: Step 1.1: Use a Butterworth filter to remove high - frequency noise signals. The Butterworth low - pass filter is represented by the following formula of the square of the amplitude with respect to frequency: Where n is the order of the filter, ω c is the cut-off frequency, ω p is the passband edge frequency. Since the vibration frequency is below 4000Hz, a filter is used to filter out high-frequency noise signals greater than 4000Hz. Step 1.2: Perform FFT on the acceleration signal {s(k), k = 1, 2, …, M} after removing high - frequency noise signals to obtain the spectrum S(i): where the length M of the spectrum = length(X), j is the imaginary unit; Step 1.3: Set the tooth passing frequency f to be removed z As shown in Equation (3), its multiple frequency is c·f z , and the initial value of c = 1 f z = Ω / 60(3) where, Ω is the spindle speed; Step 1.4: Consider the tooth passing frequency f z There are fluctuations. Introduce a floating constant a, with a default value of a = 0.

5. Set the spectral amplitude of the interval corresponding to the frequency (f z - a, f z + a) to zero, as shown in Equation (4); Step 1.5: If i is greater than c·f z +a and less than N / 2, where N is the number of points of the signal after Fourier transform, then let c = c + 1, and repeat Steps 1.3 to 1.5; otherwise, proceed to the next step; Step 1.

6. Since the spectrum after FFT is symmetric about the center, let obtain the complete spectrum after removing the tooth passing frequency and its multiples; Step 1.7: Perform inverse Fourier transform (IFFT) on the spectrum after removing the tooth passing frequency and its multiples to obtain the de - frequency signal s(k), as shown in Equation (5):

3. The multi - feature fusion method for robot milling edge chatter based on PSO - SVM - RFE according to claim 2, characterized in that, the specific steps of Step 2 are as follows: Step 2.1: Extract the time-domain features of the acceleration signal after frequency removal. Let the time series of the collected signal be {s k (t), k = 1, 2, …, M}, and the extracted time-domain feature set be T d . The features are as follows: (1) Maximum value a 1 : Describes the highest peak value of the vibration signal in the milling process, as shown in Equation (6): a 1 = max(|s k |)(6) (2) Average value a 2 : Describing the average change trend of the vibration signal in the milling process, as shown in Equation (7): (3) Peak-to-peak value a 3 : Describes the difference between the highest peak and the lowest peak of the vibration signal during the milling process, as shown in Equation (8): a 3 = a 1 - min(s k )(8) (4) Variance a 4 : Describing the severity of the amplitude fluctuation of the milling vibration signal, as shown in Equation (9): (5) Standard deviation a 5 : Describing the degree of dispersion between each signal point of the milling vibration signal, as shown in Equation (10): (6) Kurtosis a 6 : Describes the peak state of the probability distribution of vibration signal points. It is a dimensionless parameter, independent of the amplitude of the signal, but very sensitive to the impact characteristics of the signal. When a 6 = 3, it indicates that the peak value of the signal distribution is close to the peak value of the normal distribution at this time; when a 6 > 3, it indicates that the peak value of the signal distribution is greater than the peak value of the normal distribution at this time, that is, the degree of dispersion of signal points decreases, also known as positive kurtosis; when a 6 < 3, it is negative kurtosis. At this time, the degree of dispersion of the signal distribution increases, and the peak value of the signal distribution is lower than the peak value of the normal distribution curve, as shown in Equation (11): (7) Skewness a 7 : Describes the degree of asymmetry in the distribution of milling vibration signal points. Corresponding to kurtosis, kurtosis describes the situation of the peak of the signal point distribution compared to the normal normal distribution, while skewness describes the degree of symmetry of the tails on both sides of the peak of the signal distribution. Therefore, in the case of a normal distribution, a 7 = 0, that is, it is completely symmetric on both sides of the peak. When a 7 < 0, it is called left skewness, that is, at this time the signal points are mainly on the right side of the mean, showing that the left tail of the signal point distribution curve is longer; when a 7 > 0, it is called right skewness, and at this time the right tail of the distribution curve is longer. Skewness is also a dimensionless characteristic, as shown in Equation (12): (8) Root mean square a 8 : It describes the average power and average energy of the vibration signal. Different from the mean value, the root mean square value first squares and sums the data, then averages, and finally takes the square root. It can not only describe the situation of the mean value to a certain extent, but also reflect the degree of data dispersion, as shown in Equation (13): (9) Waveform factor a 9 : Describing the waveform change trend of the milling vibration signal, as shown in Equation (14): (10) Peak factor a 10 : Describing the vibration moment of the milling vibration signal, as shown in Equation (15): (11) Pulse factor a 11 : Describing the sensitivity of the milling vibration signal to the pulse, as shown in Equation (16): (12) Margin factor a 12 : Describing the impact situation received in the milling vibration signal, as shown in Equation (17): The above are the time-domain features extracted from the signal time domain, including dimensional features and dimensionless features. The initial time-domain feature set T d is as follows: T d = {a 1 , a 2 ,..., a 12} (18); Step 2.2: Extract the frequency-domain features of the acceleration signal. Different from the time-domain analysis that reflects the change trend of the signal amplitude over time, the frequency-domain analysis reflects the frequency information contained in the signal from another dimension. When chatter occurs during robot milling, not only will there be obvious changes in the signal time domain, but also chatter frequencies close to the natural frequency of the robot system will be generated in the frequency spectrum. The amplitudes of these frequencies directly indicate the instability degree of the system at this time. FFT is the main method of frequency-domain analysis, and its transformation formula is shown in Equation (2). After the signal is transformed by FFT, the frequency vector is {f i , i = 1, 2, …, M}, where f i represents the i-th frequency value of the frequency vector, and the amplitude of the frequency power spectrum corresponding to this frequency is b i , that is, the abscissa and ordinate in the frequency spectrum diagram. Let the frequency-domain feature set composed of frequency feature quantities be T f , which is specifically shown as follows: (1) Center of gravity frequency a 13 : Describes the center of gravity of the frequency in the spectrum of the robot milling signal, that is, it reflects the frequency distribution in the spectrum at a certain moment. The formula for the center of gravity frequency is shown in Equation (19): (2) Mean square frequency a 14 : Describes the change of the main frequency band of the power spectrum in the robot milling signal spectrum, as shown in Equation (20): (3) Frequency variance a 15 : Describing the discrete situation of the energy distribution of the power spectrum of the robot milling signal, as shown in Equation (21): Frequency domain analysis transforms continuous time-domain signals into discrete frequency domains through FFT, which can more accurately analyze the content contained in the signals. The frequency domain feature set T f is as follows: T f = {a 13 , a 14 , a 15} (22) Step 2.3: Perform wavelet packet decomposition on the de - frequency acceleration signal: Use wavelet packet decomposition (WPD) to decompose the signal into multiple sub - signals. At the same time, use energy entropy as the time - frequency - domain feature. When using WPD to perform multi - layer decomposition on the signal, the original signal is first decomposed by the filter into an approximate signal and a detail signal, that is, the signal passes through a low - pass filter to obtain a low - frequency signal and through a high - pass filter to obtain a high - frequency signal. Then, at the next layer, the approximate signal and the detail signal are further decomposed, and the signal is divided into two parts again, low - frequency and high - frequency. WPD is shown in Equation (23), and the signal reconstruction of each node is shown in Equation (24): Wherein, a is the high-pass filter coefficient, b is the low-pass filter coefficient, d is the wavelet packet decomposition coefficient, n is the decomposition level, j is the frequency band node number, l = k + 1, and k = M / 2 n , representing the number of wavelet packet coefficients of each frequency band in the nth layer, and M is the number of sampling points of the time series S; Step 2.4, Flutter-sensitive frequency band selection: To improve the recognition ability of the energy entropy of time-frequency domain feature quantities for the initial flutter frequency, multiple sub-signals after WPD are screened to remove the sub-signal frequency bands that are weakly sensitive to the flutter state, and the correlation coefficient ρ and the variance of the fluctuation of the frequency band energy ratio σ i (t) 2 Select the sub-signals with strong sensitivity to robot milling flutter after WPD: (1) Correlation coefficient ρ After performing FEA and denoising on the collected original acceleration signal, the main component contained in the frequency-removed signal is the flutter signal. The correlation coefficient ρ is used to reflect the correlation degree between each sub-signal s i and the original signal S, and its value range is [-1, 1], as shown in Equation (25): where s i is the i-th sub-signal, S is the original signal, cov(s i , S) is the covariance between s i and S, and are the means of s i and S respectively, D(s i ) and D(S) are the variances of s i and S respectively, and n is the number of wavelet packet decomposition levels. When |ρ(s i , S)| ≤ 0.1, it indicates that the sub-signal is slightly correlated with the original signal; 0.1 < |ρ(s i , S)| ≤ 0.5 is real correlation; 0.5 < |ρ(s i , S)| ≤ 1 is significant correlation. For the acceleration signal of robot milling, the larger the absolute value of the correlation coefficient between the reconstructed signal of the sub-signal wavelet packet coefficient and the original signal, the greater the degree of correlation and the more useful signal components it contains; on the contrary, it indicates that the signals contained in the frequency band are mainly interference noises; (2) Fluctuation variance of the frequency band energy ratio In the stable milling state of the robot, the fluctuation of the energy ratio of each sub-signal band within each sampling period T follows a normal distribution, and the energy ratio of the sub-signal fluctuates up and down around the mean. However, when chatter occurs during robot milling, the energy of some sub-signals will change sharply and finally gather at the natural frequency of the robot. The energy ratio of the sub-signal containing the chatter frequency increases, and the variance of the reconstructed signal of the wavelet packet coefficient of the sub-signal increases significantly numerically. Therefore, the variance of the energy ratio fluctuation σ i (t) 2 is used to represent the sensitivity of the reconstructed signal of the sub-signal to the abnormal chatter state, as shown in Equation (27): where μ i is the mean of the energy ratio of the i-th sub-signal, and P i (t) is the energy ratio of the i-th sub-signal at time t. The larger the variance σ i (t) 2 , the better the sensitivity of the sub-signal to the milling chatter of the robot; conversely, it indicates a weakly sensitive sub-signal to chatter. After calculating the variances of the energy ratios of each sub-signal, they are rearranged in descending order as {σ 2} = {σ 1 2 , σ 2 2 , …, σ j 2 , where j = 1, 2, …, 2 n}; (3) Definition of the chatter - sensitive frequency band According to the correlation degree between each sub - signal and the original signal and the fluctuation of the sub - signal energy ratio with the chatter state, remove the redundant interference of the weakly sensitive frequency band sub - signals to the system calculation. It is set that the chatter - sensitive sub - signals need to meet the following two conditions simultaneously: 1) The correlation coefficient is greater than the set value [ρ] ρ > [ρ] (28) 2) The fluctuation variance of the band energy ratio is among the top N i values {σ 2} = {σ 1 2 , σ 2 2 , …, σ j 2 , j = 1, 2, ..., N i}, N i ≤2 n (29) N i = round(λ·N) (30) where round is the rounding function in MATLAB and λ is the sensitivity coefficient; Step 2.5: Extract the energy entropy of the time-frequency domain feature quantity: Calculate the energy entropy for the selected sensitive frequency band. The calculation formulas for the energy entropy are shown in Equations (31) and (32): where E i is the energy of the i-th frequency band, and P i is the proportion of the energy of the i-th sub-signal after normalization to the total energy, and there is According to formulas (28) and (29), the correlation coefficient is set to [ρ]=0.1, the sensitivity coefficient λ = 0.7, N i is 11. The signal is decomposed into 16 sub-signals by using the "db3" wavelet basis for four-layer wavelet packet decomposition. According to the sensitivity coefficient and the correlation coefficient, finally, the sub-signals s2, s9, s10, s11, and s12 are removed and not considered when calculating the time-frequency domain characteristic energy entropy.

4. The multi-feature fusion method for robot milling edge chatter based on PSO-SVM-RFE according to Claim 3, characterized in that the specific steps of Step 3 are as follows: Step 3.1: Establish a support vector machine (SVM) classification model. Assume that the data set {(x i , y i ), i = 1, 2, …, n} is a set of training data, where i is the number of samples, x i is an n-dimensional feature vector, y i is the class label with a value range of {-1, 1}. The distance from the data point isolation margin L1 or L2 to the decision surface is defined as 1 / ||w||. Since the two boundaries are symmetric about the decision surface, the margin interval between L1 and L2 is 2 / ||w||. Therefore, the problem of exploring the optimal separation plane is transformed into solving the minimum ||w||. The SVM decision function is shown in formula (33): f(x i ) = ω·x i + b(33) where ω = {ω 1 ; ω 2 ; …; ω d} is the weight vector that determines the direction of the classification hyperplane. In the vector, d represents the number of eigenvalue, b is the displacement term, and the distance between the classification hyperplane and the origin is b / ||w||; The position of the classification hyperplane depends on the weight vector ω and the displacement b. Considering that the data samples are not completely linearly separable, slack variables ξ are introduced i , and the optimal classification hyperplane is obtained by solving the following optimization constraint equations where ξ i = max(0, 1 - y i (ω·x i + b)), which allows the samples to not satisfy the constraints to a certain extent. Each data sample has a corresponding slack variable, indicating the degree to which the sample does not satisfy the constraints. C is the penalty function, whose value is greater than zero, representing the trade-off between the training error and the margin. The larger the value of C, the greater the penalty for classification, the higher the fitting degree of the training samples, but the generalization ability of the model will decrease, is the mapping function; Solve Equation (34) using the Lagrange function, and the optimized constraint equation is written as follows: Among them, α i is called the Lagrange multiplier, which is a constant determined during the optimization process. To avoid an overly high dimensionality of the feature space and reduce the complexity of the classification model, a kernel function is introduced, which satisfies the Mercer theory and enables the training model to have a minimization algorithm; The SVM non-linear decision function is described by Equation (36): Among them, is the optimal vector obtained from Equation (35), In actual model training, the radial basis function is selected as the kernel function, as shown in Equation (37): In the formula, σ is the width parameter of the Gaussian kernel function, which controls the radial action range of the kernel function. The generalization ability of the training model first increases and then decreases with the increase of σ. The iterative method is used to find the appropriate parameter σ to make the training model reach the optimal generalization ability; Step 3.2: Establish a support vector machine recursive feature elimination (SVM-RFE) model for feature vector sorting: Use the cost function as the weight of each feature of the SVM. For the linear SVM, the cost function for eliminating the i-th feature is shown in Equation (38): DJ(i) = (ω i ) 2 (38) For the non-linear SVM, the cost function for eliminating the i-th feature is shown in Equation (40): where DJ(i) is the cost function of the i-th removed feature, and H is the element y i y j K(x i ,x j ) matrix, H(-i) is the matrix after eliminating the i-th feature, K is the kernel function, which represents the similarity between x i ,x j ; Step 3.3: Select the sensitive features of robot milling edge chatter based on the SVM-RFE model: Use the feature vector set composed of the 12 time-domain features, 3 frequency-domain features, and 1 time-frequency domain feature extracted in Step 2 as the input of the SVM-RFE. Before performing the SVM-RFE algorithm, first normalize the data to improve the calculation accuracy and time of the model. Use min-max normalization, as shown in Equation (41): Then, the cross-validation algorithm (K-CV) is used to select the best penalty function and kernel function parameters. The penalty function C = 3.0314, and the Gaussian kernel function g = 4. After the SVM-RFE algorithm, the feature ranking vector r = [a 16 , a 9 , a 13 , a 14 , a 15 , a 6 , a 4 , a 8 , a 5 , a 3 , a 10 , a 7 , a 12 , a 2 , a 11 , a 1 .

5. The multi-feature fusion method for robot milling edge chatter based on PSO-SVM-RFE according to Claim 4, characterized in that the specific steps of Step 4 are as follows: Step 4.1: Select the feature quantities for fusion: Construct 16 feature subsets in sequence according to the feature sorting vector obtained in Step 3. Each subset has one more feature quantity than the previous subset. Substitute the 16 feature subsets into the SVM for classification training respectively, and select the sensitive feature set [a16, a9, a13]; Step 4.2, perform multi - feature fusion based on the Particle Swarm Optimization (PSO): Before iteration, the PSO algorithm initializes a group of random particles as the solutions of the function, and then continuously iterates in the D - dimensional space to optimize the optimal position, that is, to search for the global extreme value P g , and the particles continuously update their own velocities V and positions X during the iteration process. The "flight" change formula is as follows: where ε is the inertia weight; i = 1, 2, ..., n; k is the current iteration number; V in is the particle velocity; c 1 , c 2 is the acceleration factor, and its value range is [0, 2]. c 1 determines the optimal step length of the particle itself, and c 2 determines the optimal step length of the particle globally; r 1 , r 2 is a random number, and its value range is [0, 1]; Use the PSO algorithm to fuse the sensitive feature sets [a 16 , a 9 , a 13 screened in step 4.1 into a fusion feature a17 that is strongly correlated with the flutter state. The fusion objective formula is shown in Equation (44): where a 17 is the feature a 16 , a 9 and a 13 are new features fused with certain coefficients, and the sum of the coefficients l 1 , l 2 , l 3 is 1, and the optimization objective is to maximize the cost function of the fused feature; Before optimizing the PSO-SVM-RFE algorithm, initialize c 1 = c 2 = 1.49445, k max = 100, n = 20. After optimization, the fusion coefficients l 1 = 0.8393, l 2 = 0.1566, l 3 = 0.0041.

Images