Locating Faulty Rolling Element Bearing Signal by Simulated Annealing

agerasiaetherealΤεχνίτη Νοημοσύνη και Ρομποτική

24 Νοε 2013 (πριν από 3 χρόνια και 8 μήνες)

86 εμφανίσεις


Locating
Faulty
Rolling Element Bearing
Signal by Simulated Annealing


Mid Year Report

AMSC 663,
2012


Jing Tian


Course Advisor:
Professor Radu Balan
, Professor Kayo Ide

Research
Advisor: Dr. Carlos Morillo

University of Maryland, College Park


Abstract

Vibration

acceleration

signal

is
widely
used in
the health

monitoring of rolling element
bearing
s
.

A critical work of the bearing fault diagnosis is

l
ocating the
optimum
frequency
band
that

contains faulty bearing signal
, which is usually masked by noise
.
This project diagnoses the
bearing fault with a combination of methods. The optimum frequency band is located by
maximizing spectral kurtosis with simulated annealing. Then the fau
lt feature frequency
component is extracted from the frequency band by using envelope analysis.

Modules of the

program

were

validated by a combination of analytic

work and simulation work. The
whole
program will be validated by

an open database
of

both nor
mal and faulty bearing vibration
signals.



I.

Introduction

Bearing provides relative rotational freedom and transmits a load between two structures.

It
is widely used in electromechanical systems.

Rolling element bearing
is a major source of failure

in electromechanical systems. For
example,
bearing faults account for more than 40% of th
e induction motor’s failure [1], and
g
earbox bearing failure is the top contributor of the wind turbines downtime

[2, 3]. Bearings are
inexpensive devices, but the fa
ilure of bearing is costly.
A $5,000 wind turbine bearing
replacement can easily turn into a $250,000 project, not to mention the cost of downtime [4]
.
In
1987, LOT Polish Airlines Flight 5055 Il
-
62M crashed because of failed bearings in one engine,
killin
g all
the183 people on the plane [5]. In
-
situ health monitoring is used to improve the
condition
-
based maintenance, which reduces the frequency and the loss of the bearing failure.

In the bearing health monitoring, e
arly detection of the bearing fault is a

major concern for
the industry
.
Vibration acceleration signal is widely
used in this purpose because it
is sensitive to
the bearing fault

and it can
be monitored in
-
situ
.

The objective of
the vibration signal
bearing fault detection is to 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

x(t)

is the monitored vibration signal;
s(t)
is the faulty bearing signal;
v(t)
is the noise,
which is unknown
.

An industrial practice to

test the existence of
s(t)

is to test if a unique frequency component
of
s(t)
-

the fault feature frequency component
can be extracted from
x(t)

or not.

If the fault
feature frequency component is ex
tracted, the hypothesis that the bearing is faulty is true,
otherwise the hypothesis is false.

According to the research in [6], f
aulty bearing signal

s(t)

is a modulated signal


s(t) = d(t)c(t)


where
d(t)
is the modulating signal.

It is a result of the periodic impact between the bearing’s
rolling elements and the fault on the bearing’s contact surface.

Its frequency co
mponent is the
fault feature frequency, which is illustrated in a simulated faulty bearing signal in Fig.1.
The
fr
equency is provi
ded by the bearing manufacturer;
c(t)
is the carrier signal, which is a result of
the loading and vibration transfer function. This signal is usually unknown.



Fig. 1, Faulty bearing signal
s(t)

f
Fault

is the fault feature frequency


Met
hods like envelope analysis have been developed to extract the fault feature frequency.
The problem is that i
n the presence of noise the extraction may fail.

The solution is to
band
-
pass
filter the vibratio
n signal in the frequency domain, as shown in

a si
mulated vibration signal in

Fig. 2.



Fig. 2, Vibration signal in the frequency domain


The challenge to design the filter is that the optimum frequency band to band
-
pass filter the
faulty bearing signal is usually unknown. This project provides a
solution to find the optimum
frequency band.


II.

Methodology

0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0.09
0.1
-0.5
0
0.5
1
Time(s)
Amplitude
0
2000
4000
6000
8000
10000
12000
14000
0
2
4
x 10
-4
Frequency (Hz)
Amplitude

1/
f
Fault
0
500
1000
1500
2000
2500
0
0.005
0.01
Frequency(Hz)
Amplitude
Faulty bearing
signal
This project
locates

the optimum frequency band by optimizing the band
-
pass filter with

simulated annealing

(SA).

The idea is,
t
he frequency band
dominated

by
the faulty bearing signal

is non
-
Gaussi
an, and
therefore it

has
a high

spectral kurtosis value

[7].

In the presence of white Gaussian noise, by
maximizing the SK
, the optimum frequency band for the faulty bearing signal
can be

found. The
optimization problem is to maximize SK in terms of the ce
ntral frequency, bandwidth, and the
order of the finite impulse response (FIR) band
-
pass filter.




where

f
c

is the frequency band’s central frequency; Δ
f
is the width of the band; M is the order of
FIR filter;
f
Faul
is the fault feature frequency;
f
s
is

the sampling rate.

When the optimum frequency band is obtained, envelope analysis is applied to the filtered
signal to extract the bearing faulty feature frequency.

Fig. 3 shows the flow chart of the algorithm.



Fig. 3, Flow chart of the algorithm


III.

Implementation and Validation


(1)

Filter
-
bank

x(n)

is the sampled version of the vibration signal
x(t)
. It has N points.
At first, the vibration
signal x(n) is band
-
pass filtered by a FIR filter
h

to produce the filtered signal

y(n)
:

2
2
;
2
)
,
,
(
f
f
f
f
f
f
f
to
Subject
M
f
f
SK
Maximize
s
c
s
Fault
c









FIR filter
h
i
(
f
ci
,
Δ
f
i
, M
i
)
SK
SA
Maximize SK by
f
c
,
Δ
f
, M
Optimized
FIR filter
h(
f
co
,
Δ
f
o
, M
o
)
EA
x(n)
y
i
(n)
SK
i
x(n)
y
o
(n)
FFT
a(n)
Magnitude
A(
f
)
|A(
f
)|
Maximized SK
SK
o
f=
f
Fault
?
The bearing
is faulty
The bearing
is normal
Yes
No
x(n)
is the sampled vibration signal;
y
i
(n) is filtered output of the
ith
FIR filter
h
i
;
SK
i
is the SK of the
y
i
(n);
y
o
(n) is the output of the optimized FIR filter;
a(n)
is the envelope of
y
o
(n) ;
A(
f
)
is the FFT of
a(n)




h
d
(n)
is the
impulse response of the filter



w(n) is the window function
. In this project, Hamming window will be used:



This algorithm is implemented by Matlab’s build
-
in function “fir1”.

Before optimizing the filter, initial input is obtained by calculating SK fo
r the signal filtered
by an FIR filter
-
bank.

The filter
-
bank has a
structure

of binary tree as
shown

in Fig.
4
.
S
k,j

denotes the
j
th filter at the
k
th level.

When the signal
is

processed by the filter
-
bank,
sub
-
signals
corresponding

to

the filters are

obta
ined.



Fig.
4

Structure of the FIR filter
-
bank


To validate the program, a test signal having four frequency components is used. The
sampling rate is 1024Hz, and 10240 data points are used.

The signal is shown in Fig.
5
.


)
3
/
510
2
cos(
4
)
2
/
384
2
cos(
6
)
4
/
200
2
cos(
8
)
2
2
cos(
10














t
t
t
t
x



h
n
x
n
y


)
(
)
(
)
(
)
(
n
w
n
h
h
d

)
2
/
(
]
2
/
2
/
)
2
/
sin[(
]
2
/
2
/
)
2
/
sin[(
)
(
M
n
f
f
f
M
n
f
f
f
M
n
n
h
s
c
s
c
d










M
n
M
n
n
w




0
),
2
cos(
46
.
0
54
.
0
)
(

Level 0
Level 1
Level 2
Level 3

Level k
S
k,j

0
f
s
/2

Frequency
S
0,1


S
2,3
S
2,2
S
2,1
S
1,2
S
1,1

S
3,1
S
3,2
S
3,3
S
3,4
S
3,5
S
3,6
S
3,7
S
3,8
S
2,4
Fig.
5
. Test Signal (Left: time domain; right: frequency domain; dotted red line: edge of a
frequency band)


A 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 frequency bands
. Results in the frequency domain
are illustrated in Fig. 6
.
All the frequency components locate on the right place of the frequency
axis. But the magnitude was changed.


Fig.
6
.

Frequency components obtained by applying the filter
-
bank

(Above left: pass
-
filtered
in
[0, 128]Hz
; above right: pass
-
filtered in
[
128
,
256
]Hz
;

below left: pass
-
filtered in
[
256
,
384
]Hz
;

below right: pass
-
filtered in
[
384
,
512
]Hz
)


Response of the magnitu
de is obtained by calculating the magnitude of the Fourier transform
of the filters


coefficients, as illustrated in Fig.

7.

W
e can observe that the
magnitude response

at
384Hz
is smaller than one
. That

explains why the
384Hz frequency component has
smaller
magnitude than its real value.


Fig. 7.

Magnitude response of the filter
-
bank


(2)

Spectral Kurtosis

Then the spectral kurtosis of the filtered signal
y(n)

is calculated. Spectral kurtosis is d
efined
as follows:




where
κ
r

is the
r
th order cumulant.
Y(m)
is the DFT of the signal
y(n)
:




Both y(n) and Y(m) are N points sequences. SK is a real number.

To estimate SK, the formula for joint cumulant is used:




According to [8],
DFT of a stationary signal is a circular complex random

variable, and
E[Y(m)
2
]=0, E[Y
*
(m)
2
]=0.
Therefore, we have




Program of SK is validated by a white Gaussi
an noise and a cosine function
, which are
plotted in Fig. 8.
The
sampling rate is
fs
=
1024Hz, and
the number of the

data points
is N=
1024.





Fig. 8. White Gaussian noise(Above left: time domain; above right: frequency domain), and
cosine function (Below left: time domain; below right: frequency domain)


2
*
2
*
*
4
)}]
(
),
(
{
[
)}
(
),
(
),
(
),
(
{
m
Y
m
Y
m
Y
m
Y
m
Y
m
Y
SK



1
,...,
1
,
0
,
)
(
)
(
1
0
2







N
m
e
n
y
m
Y
N
n
N
n
m
i

)]
(
)
(
[
2
)]
(
)
(
[
)]
(
)
(
[
)]
(
)
(
)
(
)
(
[
*
*
*
*
*
m
Y
m
Y
E
m
Y
m
Y
E
m
Y
m
Y
E
m
Y
m
Y
m
Y
m
Y
E
SK




2
}]
|
)
(
{|
[
}
|
)
(
{|
}]
|
)
(
{|
[
}]
|
)
(
{|
[
2
}
|
)
(
{|
2
2
4
2
2
2
2
4




m
Y
E
m
Y
E
m
Y
E
m
Y
E
m
Y
E
SK
If we investigate two signals
analytically,
we can get the following results.
Fourier transf
orm
of the white Gaussian noise also follows Gaussian
distribution
, and therefore its spectral kurtosis
is 0. For a cosine function
y=cos(
2
π
f
t)
,
f=20,

its Fourier transform is


)]
(
)
(
[
2
1
)
(
f
m
f
m
m
Y








where
δ
(m) is the delta function
.

W
hen
m=f

and
m=
-
f
,
|Y(m)|=1/2
, otherwise
|Y(m)|=0
.

SK of the cosine function is estimated with a given number of data points N=1024 as follows:

510
2
/
)
5
.
0
5
.
0
[(
/
)
5
.
0
5
.
0
(
2
]
/
)
|
)
(
|
[(
/
)
|
)
(
|
(
2
}]
|
)
(
{|
[
}
|
)
(
{|
2
2
2
2
4
4
2
1
2
1
4
2
2
4













N
N
N
m
Y
N
m
Y
m
Y
E
m
Y
E
SK
N
k
k
N
k
k


The numerical result is obtained as follows:
SK for the white noise is
-
0.0406
; SK for the
cosine function
is 510
.0000
.

Both
values

match the analytical results.



(3)

Simulated annealing

T
h
e process of estimating SK as a function of the FIR filter is optimized by simulated
annealing

(SA) [9]
, which is

a
metaheuristic global optimization tool.

The flowchart of
implementing is illustrated in Fig.
6
. In reach iteration, there is a chance that a worse case would
be accepted and thus simulated annealing can avoid the searching being trapped in a local
extremum.


The algorithm is validated by a one
dimension

function and a three dimension function.
The one dimension function tested the algorithm

s ability to locate the global optimum when the
function has many local extremum.
The function is plotted in Fig. 9. The global minim
um is y=
-
100 when x
=
π
.


2
)
(
)
cos(
100





x
x
y


Fig. 9. One dimension function



Setting of the simulated annealing is illustrated in Fig.
10
. Initial temperature T=1000.
Every time when a new result is accepted, T is reduced as T=0.99T. The initial input is a random
number between [
-
10, 10]. In every iteration a step is added to the input. The step is a random
-10
0
10
20
-100
0
100
200
x
y
number, the step size is reduced
as S=0.9999S when T drops.

The algorithm stops at 1
,
000
iterations.



Fig.
10
,
Setting of the simulated annealing for one dimension function



The minimum found was y=

-
99.999
7. The corresponding
minimize
r is x=

3.1440
.

The
searching process is illustrate
d in Fig. 11.




Fig. 11. Searching process (Left: variable; right: function value)


The algorithm was further validated by a three dimension function. This function can be
regarded as the addition of three functions

(Fig. 12)
. Global minimum is achieved

when all the
three functions


global minimum are reached. The global minimum is y=
-
150, when
x
1
= π, x
2
=15,
and x
3
=
-

4π.

)
4
/(
)
4
sin(
50
)
15
(
)
(
)
cos(
100
3
3
2
2
2
1
1












x
x
x
x
x
y


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

Fig. 12. Three dimension function


Fig.
13.

describes the setting of the algorithm.

Initial temperature T=1000. Every time when a
new result is accepted, T is reduced as T=0.99T. The initial input is a three
-
element vector. Each
element is a random number between [
-
10, 10]. In every iteration a step vector is added to the
input. Elements

of the step vector are random numbers. Amplitude of the vector is reduced as
S=0.9999S when T drops. The algorithm stops at 100,000 iterations.




Fig. 13. Setting of the simulated annealing for three dimension function



The minimum found was y=
-
149.7608
. The corresponding
minimize
r is x=[
3.179
6,
14.808
8,
-
12.60
60].

The searching process is illustrated in Fig. 14.


-20
-10
0
10
20
30
-100
0
100
200
x
y
2
1
1
)
(
)
cos(
100





x
x
y
2
2
)
15
(


x
y
)
4
/(
)
4
sin(
50
3
3






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

Fig. 14. Searching process (Left: variable; right: function value)


(4)

Envelope analysis

When the optimized frequency band is found,
envelope analysis is applied to the filtered
signal.
Th
e enveloped signal is obtained
from the m
agnitude of the analytic signal

which is
constructed
via Hilbert transform
:






Analytic signal



The envelope is the magnitude of the analytic signal



Fig.
15

shows the effect of envelope analysis on a modulated signal


Fig.
15
, Effect of envelope analysis


The program is validated by a modulated signal

in Fig. 16
:


2
1
)
1
(
x
x
y



where

)
10
2
cos(
1
t
x



)
4
/
200
2
cos(
100
2




t
x





d
t
h
y
t
y
o
o
)
(
)
(
)
(
ˆ






t
t
h

1
)
(

)
(
ˆ
)
(
)
(
t
y
j
t
y
t
y
o
o
a


|
)
(
|
)
(
t
y
t
a
a

0.032
0.033
0.034
0.035
0.036
0.037
0.038
-0.5
0
0.5
Time(s)
Amplitude

Original signal
Hilbert transform of
the original signal
Enveloped signal

Fig.

16.

Modulated signal (Left: time domain; right: frequency domain)


Frequency component of the modulating signal is expected to be extracted after the envelope
analysis and Fourier transform.

After envelope analysis,
de
-
modulated signal is obtained

(Fig. 17)
.
The de
-
modulated signal
has the same frequency as the modulating signal, and the same amplitude as the original
modulated signal.



Fig.

17.

De
-
modulated signal (Left: time domain; right: frequency domain)


IV.

Summary

All the modules of the program achieved their designated function. Some
topics can be
investigated to improve the method.
Firstly, the reduced magnitude after filtering needs to be
compensated; secondly
, the performance of the envelope analysis module in c
ase
of noise has
not been examined.

In 2013 term, all the modules will be assembled to achieve the task of bearing fault detection.
And parallel computing will be implemented to improve the program.


2012



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


2013



February

-

Complete validation



March

-

Adapt the code for parallel computing



April

-

Validate the
parallel version



May

-

Final report and presentation



References

[1] L. M. Popa, B.
-
B. Jensen, E. Ritchie, and I. Boldea, “Condition monitoring of wind
generators,” in
Proc. IAS Annu. Meeting
, vol. 3, 2003, pp. 1839
-
1846.

[2] Wind Stats Newsletter, 2003

200
9, vol. 16, no. 1 to vol. 22, no. 4, Haymarket Business
Media, London, UK

[3] 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


[4] C. Hatch, “Improved wind turbine condition monitoring using acceleration enveloping,”
Orbit, pp. 58
-
61, 2004.

[5] Plane crash information

http://www.planecrashinfo.com/1987/1987
-
26.htm


[6]
P. D. Mcfadden, and J. D. Smith, “Model for the vibration produced by a single

point defect in a rolling element bearing,” Journal of Sound and Vibration, vol. 96,

pp. 69
-
82, 1984.


[7] J. Antoni, “The spectral kurtosis: a useful tool for characterising non
-
stationary signals”,
Mechanical Systems and Signal Processing
,
20,
pp.282
-
307, 2006

[8] P. O. Amblard, M. Gaeta, J. L. Lacoume, “Statistics for complex variables and signals
-

Par
t I:
Variables”, Signal Processing
53
, pp. 1
-
13, 1996

[9] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi, "Optimization by Simulated Annealing".
Science

220

(4598), pp. 671

680, 1983