Locating Faulty Rolling Element Bearing Signal by Simulated Annealing

1. AMSC 663 Mid Year Presentation, Fall 2012 Locating Faulty Rolling Element Bearing Signal by Simulated Annealing Jing Tian Course Advisor: Dr. Balan, Dr. Ide Research Advisor: Dr. Morillo

2. Background • Rolling element bearings are used in rotating machines in different industry sections. Bearings Bearing http://en.wikipedia.org/wiki/File:J85_ge_17a_turbojet_engine.jpg Computer cooling fan Gas turbine engine Bearings inside http://en.wikipedia.org/wiki/File:Silniki_by_Zureks.jpg http://en.wikipedia.org/wiki/File:Scout_moor_gearbox,_rotor_shaft_and_brake_assembly.jpg Bearing Induction motor Wind turbine gearbox

3. Health Monitoring of Bearing • Bearing failure is a concern is a concern for many industrial sections • Bearing fault is a main source of system failure, e.g.: Gearbox bearing failure is the top contributor of the wind turbine’s downtime [1, 2]. • The failure of bearing can result in critical lost, e.g.: Polish Airlines Flight 5055 Il-62M crashed because of bearing failure [3]. http://en.wikipedia.org/wiki/File:LOT_Ilyushin_Il-62M_Rees.jpg http://en.wikipedia.org/wiki/File:DanishWindTurbines.jpg Offshore wind turbines LOT Polish Airlines Il-62M • Vibration signal is widely used in the health monitoring of bearing • It is sensitive to the bearing fault. The fault can be detected at an early stage. • It can be monitored in-situ. • It is inexpensive to acquire.

4. Project Objectives • How to detect the bearing fault?Test if the vibration signal x(t) contains the faulty bearing signal s(t) • Faulty bearing: x(t) = s(t) + ν(t) • Normal bearing: x(t) = ν(t), where v(t) is the noise, which is unknown • How to test the existence of faulty bearing signal s(t)? Check if unique frequency component of s(t) can be extracted. • Faulty bearing signal is a modulated signal : s(t) = d(t)c(t) • d(t) is the modulating signal. Its frequency component is the fault signature. The frequency is provided by the bearing manufacturer. • c(t) is the carrier signal, which is unknown. • Objective of the project: given vibration signal x(t), test if the frequency component of d(t) can be extracted.

5. Approach • Use FIR filter-bank to decompose the test signal into sub-signals. • Use spectral kurtosis (SK) to locate the sub-signal which contains faulty bearing signal s(t). SK is the kurtosis of the discrete Fourier transform of the vibration signal. The sub-signal containing faulty bearing signal has higher SK [4]. • Apply simulated annealing (SA) to optimize the frequency band of the located sub-signal. • The optimum frequency band is determined by the optimum filter. • The filter is optimized by solving the following problem: • fc is the frequency band’s central frequency; Δfis the width of the band; M is the order of FIR filter; fFaulis the fault feature frequency; fsis the sampling rate. • Perform envelope analysis (EA) to the sub-signal which has the optimum frequency band to extract the fault feature frequency (modulating frequency)

6. Flow Chart of the Algorithm SA • Maximize SK by fc, Δf, M Maximized SK • SKo x(n) yi(n) SKi FIR filter • hi (fci, Δfi, Mi) SK Optimized FIR filter • h(fco, Δfo, Mo) x(n) yo(n) a(n) A(f) EA FFT Magnitude |A(f)| • x(n) is the sampled vibration signal; • yi(n) is filtered output of the ith FIR filter hi; • SKi is the SK of the yi(n); • yo(n) is the output of the optimized FIR filter; • a(n)is the envelope of yo(n) ; • A(f)is the FFT of a(n) No f=fFault? The bearing is normal Yes The bearing is faulty

7. Filter-bank

8. Band-Pass Filter the Vibration Signal • FIR filter is defined as • M is the order of the filter. bi is the coefficient of the filter. • Matlab’s built-in function “fir1” is used. It uses the following method • hd(n) is the impulse response of the filter • w(n) is the Hamming window • Therefore, the filtered signal y(n) is a function of fc, Δf, M.

9. Initial Input for Simulated Annealing • Initial input w(fc1, Δf1, M1) is obtained by calculating SK for a binary tree of FIR filter-bank. • Output of the jth filter at level k is • Structure of the filter-bank • Frequency bands are ranked according to their SK value from large to small. Top h frequency bands are selected as the initial input.

10. Validation of the Filter-bank • The test signal has four frequency components. • One component locates at the edge of two filters in frequency domain. • Two components locates near the edge of the filter.

11. Validation of the Filter-bank • Two level filter-bank is used, which has 4 filters. For sampling rate = 1024Hz, the filters are [0, 128]Hz, [128, 256]Hz, [256,384]Hz, [384, 512]Hz • Frequency components are observed in the right frequency bands. But the magnitude is reduced. Magnitude=9.5 Magnitude=7.6 Magnitude=2.85 Magnitude=2.85, 3.8

12. Spectral Kurtosis

13. Spectral Kurtosis • Definition of spectral kurtosis where Y(m) is the DFT of the signal y(n); κr is the rth order cumulant. Both y(n) and Y(m) are N points sequences. SK is a real number. • Estimation of spectral kurtosis • DFT of a stationary signal is a circular complex random variable, and E[Y(m)2]=0, E[Y* (m)2]=0. [5]

14. Result of SK • The SK have small value for white noise and high values for periodic signal. • White noise, SK = -0.0406 • Pure periodic signal , SK = 510

15. Simulated Annealing

16. Simulated Annealing Initialize the temperature T • Simulated annealing [6] is a metaheuristic global optimization tool. • In each round of searching, there is a chance that worse result is accepted. This chance drops when the iterations increase. By doing so, the searching can avoid being trapped in a local extremum. Use the initial input vector W Compute function value SK(W) Generate a random step S Keep x unchanged, reduce T Compute function value SK(W+S) No No exp[(SK(W) -SK(W+S))/T] > rand ? SK(W+S)<SK(W) Yes Yes Replace W with W+S, reduce T Termination criteria reached? Yes End a round of searching

17. Validation of SA: 1-D Function • One dimensional function optimization • The function has a global minimum and many local minimums. • The function is used to check the fundamental of the algorithm. • The global minimum is y=-100 when x= π

18. Setting of the SA Initialize the temperature T T = 1000 W: A random number in [-10,10] Use the initial input vector W Compute function value SK(W) S: A random number in [-10,10] Generate a random step S T = 0.99T Keep x unchanged, reduce T Compute function value SK(W+S) No No exp[(SK(W) -SK(W+S))/T] > rand ? SK(W+S)<SK(W) Yes Yes Replace W with W+S, reduce T T = 0.99T Termination criteria reached? 1,000 iterations Yes End a round of searching

19. Validation of SA: Result • Result • SA ends at 1,000 iterations. The optimum is found at 943th iteration. • The minimum found is -99.9996, • Optimum variable x= 3.1389

20. Validation of SA: 3-D Function • Three dimensional function optimization • The function has three variables, it is used to validate the algorithm… • The global minimum is y=-150 when x1= π, x2=0, and x3=- π.

21. Setting of the SA Initialize the temperature T T = 1000 W: A vector, each element is a number in [-10,10] Use the initial input vector W Compute function value SK(W) S: A vector, each element is a number in [-10,10] Generate a random step S T = 0.99T Keep x unchanged, reduce T Compute function value SK(W+S) No No exp[(SK(W) -SK(W+S))/T] > rand ? SK(W+S)<SK(W) Yes Yes Replace W with W+S, reduce T T = 0.99T Termination criteria reached? 100,000 iterations Yes End a round of searching

22. Validation of SA: Setting and Result • Result • SA ends at 100,000 iterations. The optimum is found at 76,455th iteration. • The minimum found is -149.5179 • Optimum variable x= [3.1437 15.6941 -12.5659]

23. Envelope Analysis

24. Envelope Analysis • The enveloped signal is obtained from the magnitude of the analytic signal. • The analytic signal is constructed via Hilbert transform. • Hilbert transform shifts the signal by π/2 via the following formula • The analytic signal is constructed • Magnitude of the analytic signal forms the enveloped signal. Enveloped signal Original signal Hilbert transform of the original signal

25. Validation of Envelope Analysis: Test Signal • A modulated signal y is used to validate the algorithm. 200Hz 190Hz 210Hz

26. Validation of Envelope Analysis • After envelope analysis, modulating signal is obtained with a different magnitude. • After FFT, the modulating frequency component is obtained. 10Hz

27. Progress • October • Literature review; exact validation methods; code writing • November • Middle: code writing • End: Validation for envelope analysis and spectral kurtosis • December • Semester project report and presentation • February • Complete validation • March • Adapt the code for parallel computing • April • Validate the parallel version • May • Final report and presentation

28. Remaining Work at this Stage • Complete the validation of spectral kurtosis. • Put the sub-programs together to a single program. • Validate the whole program.

29. References • [1] Wind Stats Newsletter, 2003–2009, vol. 16, no. 1 to vol. 22, no. 4, Haymarket Business Media, London, UK • [2] H. Link; W. LaCava, J. van Dam, B. McNiff, S. Sheng, R. Wallen, M. McDade, S. Lambert, S. Butterfield, and F. Oyague,“Gearbox Reliability Collaborative Project Report: Findings from Phase 1 and Phase 2 Testing", NREL Report No. TP-5000-51885, 2011 • [3] Plane crash information • http://www.planecrashinfo.com/1987/1987-26.htm • [4] J. Antoni, “The spectral kurtosis: a useful tool for characterising non-stationary signals”, Mechanical Systems and Signal Processing, 20, pp.282-307, 2006 • [5] P. O. Amblard, M. Gaeta, J. L. Lacoume, “Statistics for complex variables and signals - Part I: Variables”, Signal Processing 53, pp. 1-13, 1996 • [6] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, "Optimization by Simulated Annealing". Science220 (4598), pp. 671–680, 1983