Abstract

Vibration and current signals are widely used in fault diagnosis and life prediction of electromechanical transmission systems. However, due to the complex working environment of harmonic reducer in the industrial robot, using a single signal for failure analysis or life prediction may risk to false alarming due to lacking of fault information. In the early operation stage of equipment, fault information contained in the signal is too weak to be extracted. In addition, there is certain correlation among the features, which can lead to meaningless superposition of the calculation process. Therefore, a new health index construction method is proposed integrating the current signal and vibration signal and reducing the redundancy among the features in multidomains and can effectively enhance the fault information. In addition, in view of the local optimum and slow speed caused by the random initialization of the network model, an improved life prediction method is proposed to optimize BP neural network to improve the prediction performance. The proposed method is verified by the test data of the harmonic reducer. Results show that the proposed method can predict the remaining useful life of the harmonic reducer more accurately and effectively.

1. Introduction

Industrial robots are mainly composed of reducers, controllers, servo systems, drives, and other components [1], which account for more than 70% of their production cost, and the reducers account for about 33.3%–37.8%. User’s field data show that 40%–70% of industrial robot failures are caused by reducers. The reducer, especially the harmonic reducer that directly participates in the operation at the end of the arm, is more liable to failures such as flex wheel wear and fatigue fracture. The damage of harmonic reducer will lead to decreasing of the transmission accuracy, unstable operation of the arm, and even to the failure of the industrial robot. Therefore, it is of great significance to diagnose the health status of the reducer and predict its remaining useful life (RUL).

Vibration signal is most widely used for fault diagnosis and RUL prediction of electromechanical system. Researchers have proposed many methods for preprocessing and data mining of vibration signals, such as cepstral analysis and discrete wavelet transform [26]. Zhao et al. [7] proposed a time-varying nonstationary fault feature extraction method for rolling bearings based on adaptive generalized demodulation transformation, aiming at the problem that bearing fault impact amplitude is small and easily submerged by noise at low speed. Tong and Fei [8] designed a fault diagnosis method of rolling bearings based on wavelet packets and a rough set neural network, to solve the problems of low efficiency and slow speed of fault diagnosis of locomotive rolling bearings. In order to realize the intelligent diagnosis of planetary gearboxes, Hu et al. [9] proposed an intelligent fault diagnosis method based on empirical mode decomposition and a deep convolutional neural network. Zhou and Yu [10] proposed a model based on a deep neural network that can effectively extract the key features of vibration signals and improve the performance of the recognizer.

In the electromechanical system, the failure or damage on the mechanical transmission component will cause eccentricity or torque fluctuation to the motor connected, and the magnetomotive force of the motor will generate phase modulation to compensate for the torque fluctuation. Then, the fault information can be generated in the stator current signal, reflecting the change of the external load of the motor [11]. Therefore, some researchers developed motor current signature analysis (MCSA) to study the correspondence between current characteristics and faults of mechanical components. Compared with the vibration signal, the current signal is not limited to the location of the sensor, nor does it interfere with the system of the industrial robot.

However, the current signals often have strong noise that submerges the fault information, which makes the current nonstationary, and traditional signal processing methods based on Fourier transform is no longer applicable. In this sense, researchers proposed short-time Fourier transform, wavelet transform, empirical mode decomposition, etc., to analyze and deal with such nonstationary signals [1219]. Eren and Devaney [15] analyzed the transient starting current of the induction motor through discrete wavelet transform, compared the sub-bands of the state before and after the bearing failure, and identified the resonant frequency of the bearing when the motor started to detect the bearing failure. Singhet and Kumar [16] found that the continuous wavelet transform was more easily method detecting bearing faults in mechanical equipment than the spectrum analysis, and the feasibility of MCSA technology was verified by analyzing the stator current and vibration signals. Ozgonenel and Karago [17] found that the selection of the wavelet basis determines the accuracy of the analysis results, but the method has certain difficulties in practical application. Huang et al. [18] proposed an empirical mode decomposition method that was superior to wavelet transform in waveform reconstruction and was suitable for processing nonlinear signals. Saadi et al. [19] used the entropy of the first four-order natural mode function components as the health index to detect motor failure. Although the above-given methods decompose the current signal accurately, they do not preprocess the noise interference in the current signal, which will affect the accuracy of the subsequent analysis results.

Using vibration signal or current signal alone for RUL prediction will increase the fault rate of detection and false alarm due to signal loss and signal interference. In the early service stage, the fault information contained in the signal is weak and the fault characteristics of the signal are difficult to extract. Wang and Li[20] developed a reliability model extracting same features from multisignals and verified that the reliability is more accurate using more than one signal. In view of this, this paper proposes a new health index construction method using motor current signal and vibration signal. Firstly, maximum correlated kurtosis deconvolution (MCKD) and complementary ensemble empirical mode decomposition (CEEMD) are combined to improve the signal decomposition accuracy, then copula is applied to eliminate redundant features in the intrinsic model function (IMF) component feature parameters. After that, the kernel principal component analysis (KPCA) is introduced to merge multivariate feature data. The dimension reduction method is used to obtain a few features that can indicate the degradation process of the harmonic reducer and take them as degradation index. The proposed health indicator construction method is verified in case study by the experiment data of the harmonic reducer.

2. Fault Features of Harmonic Reducer

When the outer ring, inner ring, or the rolling element of the flexible thin-walled bearing fails, the faulty part and the nonfaulty part collide periodically, and the number of collisions between the faulty part and the nonfaulty part per unit time is called the fault characteristic frequency fc. The collision causes the load torque fluctuation, which further results in the electromagnetic torque and generates a nonlinear signal in the stator current. The frequency of torque variation is also fc because they are from the same source. The motion equation of the motor is as follows:where Te is the electromagnetic torque, TL is the load torque, J is the moment of inertia, and B is the friction coefficient. When the servomotor works in a steady state, the load torque TL is a constant, and the rotor speed ωr also remains unchanged. Then, we can obtain the following equations:

If there is a fault on the reducer, then the load torque TL will have a balanced state with a certain fluctuation, and in this case, the DC and cosine components can be described as follows:where Teosc, TLosc, and ωrosc are the variations caused by the load torque fluctuation. The electromagnetic torque Te, the load torque TL, and the rotor speed ωr can be expressed as follows:where Te0, TL0, and ωr0 are the electromagnetic torque, the load torque, and the rotor speed under steady state, respectively; fc is the fault characteristic frequency; Tec, TLc, and ωrc are the amplitude of the variations, φec, φLc, and φrc are their corresponding phases. According to [21], the current can be considered to contain a torque component iq and an excitation component id, and iq can be derived using the following equation:where iq0 is the torque component current in a steady state, and iqc is the variation caused by the torque fluctuation due to the reducer faults. As a result, the fault on the reducer is reflected by the current of the servomotor connected, which is shown in Figure 1.

3. Construction of Health Indicators

A new health index construction method is proposed in this study, which mainly includes three modules: data acquisition, signal preprocessing, and degradation index construction, as shown in Figure 2.

After collecting the current signal and vibration signal, it is necessary to eliminate the high-frequency noise components in the signal through filtering processing. In this case, the median filtering method is used to process the original signal for it has a good performance for high-frequency impulse noise and can protect the edge of the signal from being blurred. Then, we use MCKD to retain the fault impact component in the signal. By taking the correlation kurtosis as the optimization object, which can accurately reflect the strength of the fault shock signal. The expression of the correlation kurtosis is as follows:where x(n) is the signal sequence, N is data points, T is the pulse period, M is the number of cycles, and m = 1,2, …, M. We also use CEEMD to smoothly decompose the nonstationary and nonlinear signals and reflect the essential characteristics of the fault in the signal. The CEEMD method mainly eliminates residual noise and improves decomposition accuracy by adding pairs of positive and negative white noise to the original signal. The specific calculation steps can be found in [22].

In addition, the parameters of the reducer may show completely different trends under the same degree of damage, which may easily lead to misjudgment of faults. The incipient defect containing in the raw signals is often too weak to be extracted, and according to [23], the degradation information of the mechanical component can be reflected in time domain, frequency domain and time-frequency domain at various levels, and different features belonging to the three domains will also contain different kinds of fault information. Therefore, the prediction result of the remaining useful life may be inaccurate or even cannot detect the incipient defect. According to [24], features from time domain, frequency domain, and time-frequency domain may not be equally sensitive to the fault and some of them can be dependent on each other. In this case, we use Copula to select the features with more fault information and reduce the redundancy. We have also summarized the feature extraction method to provide the choice of how to select features in different situations on pages 9-10 in Table 1.

In this study, we propose a health index construction method for predicting the RUL of harmonic reducers of the industrial robot. This method integrates the fault information of current signal and vibration signal. The correlation among the features is considered when constructing health index with multi-feature, and copula is applied to evaluate the correlation among the features and remove the redundant ones to reduce the redundant and improve the efficiency. The KPCA method is applied to fuse the obtained multidimensional features to obtain a health indicator that can reflect the degradation process of the harmonic reducer. The copula function is defined as follows.

The function C:[0, 1]n ⟶ [0, 1] is called the n-dimensional copula function if the following conditions are satisfied:(1)C is incremented for each variable(2)For any , satisfy: (3)For any , satisfy: (4)For all , , and , we have

Especially, for the binary copula function: for any , such that and , there is

4. RUL Prediction

BP neural network is the most widely used multilayer feedforward model; it has the advantages of nonlinear mapping and strong robustness [4246]. However, the initial weights and thresholds of BP algorithm are obtained by random initialization, and during the learning process, the weight and the threshold are dynamically, which may cause local optimum and low efficiency. This paper proposes an improved BP neural network based life prediction method, in which beetle antennae search (BAS) algorithm is used to optimize the initial weights and thresholds before train BP network, which greatly improves network performance. The framework of the BAS in BP network is shown in Figure 3.

The advantage of the BAS algorithm is that it only needs to search for a single individual and does not need to know the specific form of the objective function to achieve efficient optimization. The specific modeling steps are as follows [47]:(1)Create a random vector of the direction for the beetle antennae and normalize it:where rands(·) is the random function, and k is the spatial dimension.(2)Create the spatial coordinates of the left and right antennae that mimic the searching behavior of the beetle’s left and right antennae:where xrt is the position of the right antennae of beetle at tth iteration, xlt is the position of the left antennae of beetle at tth iteration, xt is the position of the centroid of the beetle at tth iteration, and d0 is the distance between the two antennae.(3)According to the fitness function, the smell intensity of the left and right antennae is judged; i.e., the intensity of f(xlt) and f(xrt), and the function f(·) are the fitness function.(4)Create a iterative model that mimics the search behavior of beetles:where xt+1 is the spatial position of the beetle at (t + 1) th iteration, xt is the spatial position of the beetle at tth iteration, is the step length of tth beetle detection movement, and sign(·) is the symbolic function.

In this paper, the correlation coefficient R, root mean square error (RMSE), cumulative relative accuracy (CAR), and mean absolute error (MAE) are used to evaluate the accuracy of RUL prediction model. A smaller RMSE and MAE in the experiment means that the model has a better fitting effect, a larger CAR means a better prediction accuracy, and a larger R means that the model has a stronger adaptability. The expressions of them are as follows:where n represents the total number of monitoring points, RULreal,i is the actual life of at ith monitoring point, RULest,i represents the RUL at ith monitoring point, and RULmean is the average life of the component.

5. Case Study

The proposed methods are verified by the experimental data of three harmonic reducers including current signals and vibration signals. The details of the experiment are also specified.

5.1. Accelerated Life Test of Harmonic Reducers

As shown in Figure 4, the reducer test bench mainly includes of drive motor, harmonic reducer, torque speed sensor, acceleration sensor, and load motor. The reducer test bench can conduct the accelerated life test of the harmonic reducer under different working conditions. Accelerated testing is usually adopted to speed up failures of highly reliable products for reliability analysis. In such tests, products are exposed to harsher-than-normal conditions (e.g., load, strain, temperature, voltage, vibration, and pressure) in an effort to collect failure time data in a short amount of time without changing the failure mechanisms [48, 49]. A lot of existing research using accelerating test data to verify their prognostic methods [5052], and their studies show that features extracted in accelerating test data can be used to validate the remaining useful life, indicate the degradation process, and evaluate the reliability. In our case, we need to make sure that the test of harmonic reducers cannot exceed the tolerance limit [53]. In this case, we use the stress limit of the flexible gear of harmonic reducer to be the acceleration stress level limit of accelerated testing [54]. Our test results show that the failure modes in the three experiments were all fatigue fracture, which are in accordance with the result in [55] and the engineering practice.

In this case, we set the output torque of the load motor to different multiples of the rated torque of the reducer and set the speed of the drive motor to be 2000 rpm. The current signal in the experiment is measured by Fluke 80i-110s AC/AD current clamps and collected by NI USB-628 data acquisition card, and the vibration signal is obtained using an acceleration sensor and data acquisition card. The detailed test information of the sample data are shown in Table 2, including the corresponding working conditions, sampling interval, total number of data samples, lifetime, and failure mode of each harmonic reducer sample.

5.2. Construction of Health Indicators

We use the experimental data of sample #3 as the research object for detailed analysis. The original current signal is processed by the median filter method. The current waveform before and after filtering are shown in Figure 5, and the vibration waveform before and after filtering are shown in Figure 6. It can be seen from these figures that the filtered signal merely contains high-frequency sampling noise.

The preprocessed current signal and vibration signal are decomposed by CEEMD. The 16-dimensional features of the IMF-1 component in the time domain, frequency domain, and entropy domain are extracted and shown in Table 3. Copula is used to calculate the correlation between feature dimensions to remove redundant features.

Before applying copula, the feature which can best represent the degradation trend of the harmonic reducer is selected from the above 16 features. This feature is used as the benchmark of Copula function and compared with other features, respectively, to eliminate the features with strong correlation, i.e., high redundancy. Figure 7 shows the trend curve of each feature parameter. So RMS is selected as the reference feature for analysis. The copula is used to compare the correlation between the rest of the features and RMS. If the correlation is too large, it indicates that the feature is redundant and should be eliminated (this paper used 0.3 as the correlation threshold for analysis). Results are shown in Figure 8. Among them, the three feature quantities of P3, P4, and P8 are highly correlated with RMS (annotation: feature dimension P2 is not compared with itself, so the correlation is expressed as 0).

The KPCA method is applied for data fusion and dimensionality reduction. The optimum condition is reached when the kernel parameter c = 17. In this case, the number of pivot elements is 4, which is the minimum number of pivot elements required to reach 85% of the cumulative contribution rate. The first four pivot elements are multiplied by their respective contribution rates and sum them up as the degradation index. In order to observe the change of degradation index with the degradation degree of the harmonic reducer, the mean value of the degradation index in the whole life cycle is standardized.

5.3. RUL Prediction

There are 1242 sample points in the whole life cycle of #3 harmonic reducer, which are divided according to the ratio of training set and test set 7 : 3. A total of 870 sample points are randomly selected as the training set of the BPNN model, and the remaining 372 sample points are used for testing. For the BPNN model, the input layer is KPCA1, KPCA2, KPCA3, and KPCA4 after the fusing the two signals, and the output layer is the RUL prediction value. Parts of the data are shown in Table 4. The training parameters of the BP network are shown in Table 5. The number of the hidden layer is between 4 and 13, and the root mean square error and the correlation coefficient are used as the detection targets to obtain the training parameters of the BPNN network. When the number of hidden layers is equal to 12, the RMSE is the smallest and the correlation coefficient is the largest, so the BPNN network structure is selected as 4-12-1.

The optimized genetic algorithm base BP network model (GA-BP) and BAS-BP network model are constructed by using genetic algorithm and the beetle search algorithm to optimize the initial value and threshold of BP neural network. The randomly selected 870 sample points are used for training the BP model, the GA-BP model, the BAS-BP model, and the remaining 372 sample points are used for testing. Results are shown in Figures 911.

R represents the correlation coefficient between real value and predicted value. It can be seen from the three figures that the value of R in the training set is 0.9652 and the value of R in test set is 0.9593 in BPNN model. In the GA-BP model, the value of R in training set is 0.9854 and the value of R in test set is 0.9872. In the BAS-BP model, the value of R in training set is 0.9903 and the value of R in test set is 0.9937. We can see that the BAS algorithm greatly optimizes the fitting effect of the BPNN model.

From Figure 12, we can see that when the vibration signal is used alone, the value of R of the training set is 0.9778 and is 0.9558 of the testing set. In Figure 13, when the current signal is used alone for prediction, the value of R of training set is 0.9785 and is 0.9896 of the testing set. According to Figures 1113, the prediction accuracy of single signal is lower than that of the fused signals. By comparing the experimental results, it is demonstrated that two-signal fusion can superimpose the strength of the fault information and improve the prediction accuracy.

The evaluation results of the test set are shown in Table 6. We can see that the new index combining vibration signal and current signal can better characterize the degradation process of the reducer.

We can also analyze the experiment data of harmonic reducers nos. 1 and 2. The RUL results are shown in Figures 14 and 15.

From the BAS-BP prediction results of the three sets of harmonic reducer test samples, we can see that the proposed method for constructing health index is applicable and effective to the fault characteristics of harmonic reducers.

6. Conclusion

This study proposes a health index construction method for harmonic reducers of the industrial robot, in which MCKD and CEEMD algorithms are combined to enhance the fault information, copula function is used to remove the redundant features of the IMF component, and KPCA is used to fuse the fault information of the vibration signal and the current signal to contain enough fault information. An improved RUL prediction method is also developed to verify the proposed method. An experiment of harmonic reducers is performed and the vibration and current data are obtained. The results show that the RUL is more accurate after the fusion of the current signal and the vibration signal than that of using a single signal.

In the future, we will study the current signal features of the reducer under various fault types and study the mapping method between the accelerated life test and real life.

Acronyms

MCSA:Motor current signature analysis
RUL:Remaining useful life
MCKD:Maximum correlated kurtosis deconvolution
CEEMD:Complementary ensemble empirical mode decomposition
IMF:Intrinsic model function component
KPCA:Kernel principal component analysis
BAS:Beetle antennae search algorithm
BAS-BP:Beetle antennae search algorithm optimization BP neural network
GA:Genetic algorithm
GA-BP:Genetic algorithm optimization BP neural network
ALT:Accelerating life test
DBN:Deep belief network
RMSE:Root mean square error
CAR:Cumulative relative accuracy
MAE:Mean absolute error.

Notation

Te:The electromagnetic torque
TL:The load torque
J:The moment of inertia
B:The friction coefficient
ωr:The rotor speed
Teosc:Electromagnetic torque variations caused by load torque fluctuation
TLosc:Load torque variations caused by the load torque fluctuation
ωrosc:Rotor speed variations caused by the load torque fluctuation
Te0:Electromagnetic torque in the steady state
TL0:Load torque in the steady state
ωr0:Rotor speed in the steady state
fc:Fault characteristic frequency
Tec:Teosc’ amplitude of variations
TLc:TLosc’ amplitude of variations
ωrc:ωrosc’ amplitude of variations
φec:The corresponding phase of Teosc
φLc:The corresponding phase of TLosc
φrc:The corresponding phase of ωrosc
iq:The torque component
id:The excitation component
x(n):The signal sequence
N:The data points
T:Pulse period
M:The number of time cycles, m= 1,2, …, M
k:The spatial dimension
xrt:The position of the right antennae of beetle at tth iteration
xlt:The position of the left antennae of beetle at tth iteration
xt:The position of the centroid of beetle at tth iteration
d0:The distance between two antennae
f(·):Fitness function
x t+1:The spatial position of the beetle at (t+ 1) th iteration
:The step length of tth beetle detection movement
sign(·):Symbolic function
n:Total number of monitoring points
RULreal,I:The actual life at ith monitoring point
RULest,I:The RUL at ith monitoring point estimation
RULmean:Average life of the component.

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was supported by the National Key R&D Program of China (Grant no. 2022YFB4702100), the National Natural Science Foundation of China (Grant no. 72001069), the Hebei Provincial Natural Science Foundation (Grant no. E2021202094), the Natural Science Foundation of Tianjin (Grant no. 21JCZDJC00730), and by the Hebei Provincial Youth Talent Support Program (Grant no. BJK2023031).