# Screw Remaining Life Prediction Based on Quantum Genetic Algorithm and Support Vector Machine

Authors: Xiaochen Zhang; Dongxiang Jiang
Journal: Shock and Vibration
 (2017)

---

## Abstract

To predict the remaining life of ball screw, a screw remaining life prediction method based on quantum genetic algorithm (QGA) and support vector machine (SVM) is proposed. A screw accelerated test bench is introduced. Accelerometers are installed to monitor the performance degradation of ball screw. Combined with wavelet packet decomposition and isometric mapping (Isomap), the sensitive feature vectors are obtained and stored in database. Meanwhile, the sensitive feature vectors are randomly chosen from the database and constitute training samples and testing samples. Then the optimal kernel function parameter and penalty factor of SVM are searched with the method of QGA. Finally, the training samples are used to train optimized SVM while testing samples are adopted to test the prediction accuracy of the trained SVM so the screw remaining life prediction model can be got. The experiment results show that the screw remaining life prediction model could effectively predict screw remaining life.


---

## Body

## 1. Introduction

NC machine tool is the significant production basis or fundamental unit of other equipment in manufacturing industry; thus the reliability of NC machine tool will seriously influence or restrict the development of manufacturing industry [1–3]. Ball screw, the key component of NC machine tool, is directly related to the reliability of NC machine tool. With the continuous innovation of modern equipment management technology, concepts such as condition-based maintenance (CBM) and prognostics and health management (PHM) have been put forward in order to ensure the safe operation of the equipment [4–7]. Therefore, by monitoring the working condition and predicting remaining life of the ball screw, the optimal maintenance plan of equipment can be determined so as to improve the utilization rate and security of NC machine tool.In the process of working, the vibration signals of ball screw are acquired through sensors. Then combined with sensitive feature extraction methods and life prediction modeling technologies, ball screw remaining life prediction can be achieved [8, 9]. Taking into account the fact that the monitoring data of ball screw life cycle are difficult to obtain, it is necessary to study the remaining life prediction of ball screw under the condition of finite samples. To solve the classification, regression, or control problems [10, 11], SVM maps the samples which cannot be linearly segmented from low-dimensional space to high-dimensional space, so a unique optimal hyperplane is obtained in high-dimensional space. It is worth noting that the selection of kernel function parameter will have a great effect on SVM model results during the mapping process. At the same time, the penalty factor also directly affects the results when solving the regression problem. Therefore, the optimal selection of kernel function parameter and penalty factor should be considered when SVM is used for regression prediction [12–15]. At present, using cut-and-trial method as parameters optimization strategy is not ideal. To the optimization of SVM parameters, this paper adopts QGA, which has the advantages of high efficiency and avoiding local optimum [16, 17], to search best initialized parameters of SVM, that can improve the performance of SVM and increase the stability of remaining life predicting system.To predict the remaining life of ball screw, a kind of method based on QGA and SVM is proposed here. According to the proposed method, original features are constructed with wavelet energies of all subfrequency bands from different acceleration sensors. Furthermore, Isomap is applied to reduce the dimensions of original features so the sensitive feature vectors are obtained. Then QGA is used to search the optimal kernel function parameter and penalty factor of SVM. Moreover, based on the QGA searching results, optimized SVM are trained and tested with samples, so the screw remaining life prediction model can be got. Finally, data from screw accelerated test bench are used to inspect the effectiveness of the proposed method.


## 2. Sensitive Feature Vectors

### 2.1. Screw Accelerated Test Bench

Figure1 shows the structure of screw accelerated test bench. A Beckhoff servo motor is applied to drive the testing ball screw, while the screw nut is connected to the back surface of the sliding platform. On the upper surface of the sliding platform, a rack and gear mechanism is mounted. A magnetic powder brake, the load device, is connected with the gear through a torque sensor. By the rack and gear mechanism, the magnetic powder brake could impose an axial force on ball screw. Since the torque sensor and the magnetic powder brake form a semiclosed loop control system, the axial force of ball screw can be adjusted. Besides the torque signal, acceleration signals are also acquired in the online monitoring system. One three-way accelerometer is installed on screw nut and two one-way accelerometers are installed on two bearing chocks. Three Beckhoff EL3632 modules are used to acquire acceleration signals and connected with an industrial PC. In the accelerated test, the sampling frequency is set to 5 kHz and the vibration data are stored in the industrial PC. The testing ball screw works in horizontal installation state, while it reaches 40 mm in axis diameter and 5 mm in lead. One end fixed and the other end floated installation method is applied. The initial preload force of the testing ball screw is 5% of the rated dynamic load, while grease lubrication is adopted. At the beginning of the accelerated test, the precision of the testing ball screw is grade 3. During the accelerated test, a grating scale installed on the back surface of the sliding platform is used to measure the precision of the testing ball screw. If the precision of the testing ball screw is less than grade 5, it means that the testing ball screw is faulted and the accelerated test should be stopped. The definition of three-dimensional coordinate axes is also shown in Figure 1.Figure 1
Screw accelerated test bench.


### 2.2. Feature Extraction and Dimensionality Reduction

When the ball screw is in the process of working, the surface damage such as pitting and peeling would produce pulsed force, so periodic shock vibration will occur. Wavelet packet decomposition is equipped with good time-frequency resolution and transient detection ability, thus being suitable for processing nonstationary signals such as periodic shock vibration. During the process of decomposing a signal, the output of thejth layer can be expressed as d
l
j
,
n, while l is the time serial number. Then the signal decomposition formula can be formed as(1)
d
l
j
+
1,2
n
=
∑
k
=
-
∞
∞
h
0
2
l
-
k
d
k
j
,
n
,
d
l
j
+
1,2
n
+
1
=
∑
k
=
-
∞
∞
h
1
2
l
-
k
d
k
j
,
n
,where d
l
j
+
1,2
n and d
l
j
+
1,2
n
+
1 are the outputs of the (j
+
1)th layer.Each sensor signal is decomposed by wavelet packet decomposition, and the wavelet energy of each subband is calculated. Wavelet energies of all subfrequency bands from different acceleration sensor channels construct the original features. With the purpose of reducing the dimension of the original features, Isomap, a dimension reduction method for manifold learning, is applied to extract sensitive features. To get the sensitive feature vectors, the main stages of Isomap are as follows [18, 19]:(1)
construction of the neighborhood graphG: for example, the observed data are V
(
x
i
∈
V
). If x
j is one of the nearest neighbors of x
i, G might contain the edge x
i
x
j when |
x
i
-
x
j
|
<
ε;
(2)
computation of the Euclidean length: the shortest paths are calculated for all pairs of data points;
(3)
data embedment: a new embedment of the data in Euclidean space is searched with the method of Multidimensional Scaling (MDS).


## 2.1. Screw Accelerated Test Bench

Figure1 shows the structure of screw accelerated test bench. A Beckhoff servo motor is applied to drive the testing ball screw, while the screw nut is connected to the back surface of the sliding platform. On the upper surface of the sliding platform, a rack and gear mechanism is mounted. A magnetic powder brake, the load device, is connected with the gear through a torque sensor. By the rack and gear mechanism, the magnetic powder brake could impose an axial force on ball screw. Since the torque sensor and the magnetic powder brake form a semiclosed loop control system, the axial force of ball screw can be adjusted. Besides the torque signal, acceleration signals are also acquired in the online monitoring system. One three-way accelerometer is installed on screw nut and two one-way accelerometers are installed on two bearing chocks. Three Beckhoff EL3632 modules are used to acquire acceleration signals and connected with an industrial PC. In the accelerated test, the sampling frequency is set to 5 kHz and the vibration data are stored in the industrial PC. The testing ball screw works in horizontal installation state, while it reaches 40 mm in axis diameter and 5 mm in lead. One end fixed and the other end floated installation method is applied. The initial preload force of the testing ball screw is 5% of the rated dynamic load, while grease lubrication is adopted. At the beginning of the accelerated test, the precision of the testing ball screw is grade 3. During the accelerated test, a grating scale installed on the back surface of the sliding platform is used to measure the precision of the testing ball screw. If the precision of the testing ball screw is less than grade 5, it means that the testing ball screw is faulted and the accelerated test should be stopped. The definition of three-dimensional coordinate axes is also shown in Figure 1.Figure 1
Screw accelerated test bench.


## 2.2. Feature Extraction and Dimensionality Reduction

When the ball screw is in the process of working, the surface damage such as pitting and peeling would produce pulsed force, so periodic shock vibration will occur. Wavelet packet decomposition is equipped with good time-frequency resolution and transient detection ability, thus being suitable for processing nonstationary signals such as periodic shock vibration. During the process of decomposing a signal, the output of thejth layer can be expressed as d
l
j
,
n, while l is the time serial number. Then the signal decomposition formula can be formed as(1)
d
l
j
+
1,2
n
=
∑
k
=
-
∞
∞
h
0
2
l
-
k
d
k
j
,
n
,
d
l
j
+
1,2
n
+
1
=
∑
k
=
-
∞
∞
h
1
2
l
-
k
d
k
j
,
n
,where d
l
j
+
1,2
n and d
l
j
+
1,2
n
+
1 are the outputs of the (j
+
1)th layer.Each sensor signal is decomposed by wavelet packet decomposition, and the wavelet energy of each subband is calculated. Wavelet energies of all subfrequency bands from different acceleration sensor channels construct the original features. With the purpose of reducing the dimension of the original features, Isomap, a dimension reduction method for manifold learning, is applied to extract sensitive features. To get the sensitive feature vectors, the main stages of Isomap are as follows [18, 19]:(1)
construction of the neighborhood graphG: for example, the observed data are V
(
x
i
∈
V
). If x
j is one of the nearest neighbors of x
i, G might contain the edge x
i
x
j when |
x
i
-
x
j
|
<
ε;
(2)
computation of the Euclidean length: the shortest paths are calculated for all pairs of data points;
(3)
data embedment: a new embedment of the data in Euclidean space is searched with the method of Multidimensional Scaling (MDS).


## 3. Remaining Life Prediction

### 3.1. Support Vector Machine

In this paper, SVM, which is suitable for studying small sample problems, is taken as the theoretical basis of modeling. As a machine learning algorithm developed from statistical learning theory, SVM maps inseparable learning samples from low-dimensional space into high-dimensional space through a kernel function so as to obtain an optimal hyperplane. Figure2 shows the structure of SVM. If x
i
(
i
=
1,2
,
…
,
N
) are the training samples, y
i are the training goal associated with x
i, and then searching the optimal hyperplane is transformed into solving the convex quadratic equation which can be defined as [20–22](2)
M
i
n
i
m
i
z
e
:
1
2
w
2
Subject  to:
y
i
w
·
x
i
+
b
≥
1
,
i
=
1,2
,
…
,
N
,where w and b are the coefficients of the hyperplane equation.Figure 2
The structure of SVM.According to the solution of convex quadratic programming in the optimization theory, Lagrange function is constructed as(3)
L
w
,
α
,
b
=
1
2
w
2
-
∑
i
=
1
N
α
i
y
i
w
·
x
i
+
b
+
∑
i
=
1
N
α
i
,where α
i is the Lagrange multiplier.To be applied in regression analysis, penalty factorC is introduced in the objective function. Then (2) becomes(4)
M
i
n
i
m
i
z
e
:
∑
i
=
1
N
α
i
-
1
2
∑
i
=
1
N
∑
j
=
1
N
α
i
α
j
y
i
y
j
x
i
·
x
j
K
x
i
,
x
j
Subject  to:
∑
i
=
1
N
α
i
y
j
=
0
,where the range of α
i is [
0
,
C
], i
=
1,2
,
…
,
N.In SVM model, different kernel functions and kernel function parameters have a great impact on the calculation results during the process of mapping the low-dimensional space to high-dimensional space. Meanwhile, the value of penalty factorC also determines the regression effect of SVM. To obtain the best prediction result of SVM, it is significant to choose suitable kernel function parameter and penalty factor C for SVM. Usually, the values of the kernel function parameter and penalty factor C are based on experience. Thus it is difficult to guarantee that the prediction result is the best. Therefore, QGA, a kind of heuristic algorithm, is applied to search the optimal kernel function parameter and penalty factor so as to obtain the best regression prediction effect for SVM model.


### 3.2. Quantum Genetic Algorithm

The main procedure of SVM parameters optimization can be described as follows: according to the variation ranges of kernel function parameter and penalty factor, constructing chromosomal genes with quantum bit coding system, so the evolutionary population which includes several chromosomes can be generated. Then realizing the continuous population evolution by means of cross, variation and quantum rotation gate. Finally, the optimized kernel function parameters and penalty factor can be obtained. Figure3 shows the algorithm flowchart of SVM parameters optimized by QGA.Figure 3
Algorithm flowchart of SVM parameters optimized by QGA.


#### 3.2.1. Calculation of Fitness Function

According to the quantum bit coding method, arbitrary chromosome of the populationQ
(
t
0
) can be defined as [23, 24](5)
q
j
t
0
=
α
11
t
0
β
11
t
0
∣
α
12
t
0
β
12
t
0
∣
⋯
⋯
∣
α
1
k
1
t
0
β
1
k
1
t
0
∣
α
21
t
0
β
21
t
0
∣
α
22
t
0
β
22
t
0
∣
⋯
⋯
∣
α
2
k
2
t
0
β
2
k
2
t
0
∣
⋯
⋯
∣
α
m
1
t
0
β
m
1
t
0
∣
α
m
2
t
0
β
m
2
t
0
∣
⋯
⋯
∣
α
m
k
m
t
0
β
m
k
m
t
0
,where q
j
t
0 is thejth chromosome of the t
0 generation and m is the number of the chromosomal genes. The quantum bit number of each chromosomal gene is shown by k
1
,
k
2
,
…
,
k
m.Then population ofN chromosomes can be expressed as(6)
Q
t
0
=
q
1
t
0
,
q
2
t
0
,
…
,
q
N
t
0
.In this paper, SVM kernel function parameter and penalty factor are corresponding to the genes of chromosome q
j
t
0. Since each chromosome should express same initialization state, the initial probability amplitude α
,
β
T is 1
/
2
,
1
/
2
T. As a kind of validation method for assessing the results of a statistical analysis, cross validation is applied to construct the fitness function, so the overfitting problem can be avoided. The key steps of cross validation is to segment a sample of data into several subsets, performing the analysis on one subset which is called the training set, and validating the analysis on the other subset which is called the testing set. The schematic diagram of cross validation is shown in Figure 4.Figure 4
The schematic diagram of cross validation.Then the fitness function can be formed as(7)
f
=
-
1
s
∑
p
=
1
s
M
p
,where s is the number of iterations. M
p represents the mean squared error (MSE) of the prediction results forpth testing set.


#### 3.2.2. Quantum Cross and Variation

Quantum cross means exchanging the corresponding genes of any two chromosomes in order to form two similar chromosomes. The purpose of the quantum cross is to exchange information between two possible solutions. The schematic diagram of quantum cross is shown in Figure5.Figure 5
The schematic diagram of quantum cross.By changing one or several genes of any chromosome, a new chromosome can be formed, so as to ensure the diversity of the chromosomes. This process is called quantum variation which is a good way to change the current evolution direction and avoid local optimum.


#### 3.2.3. Quantum Rotation Gate

Quantum gate is adopted to be the actuator of the evolution process in QGA. The quantum bit probability amplitude is updated based on quantum rotating angleθ
(8)
α
′
β
′
=
c
o
s
θ
-
s
i
n
θ
s
i
n
θ
c
o
s
θ
α
β
,where α
′
β
′
T is the new quantum bit probability amplitude and α
β
T is the quantum bit probability amplitude before updating.The adjustment strategy [25] of the quantum rotating angle θ is shown in Table 1.Table 1
The adjustment strategy of the quantum rotating angle.
x
i
b
i
f
(
x
)
>
f
(
b
)
Δ
θ
i
s
(
α
i
,
β
i
)
α
i
β
i
>
0
α
i
β
i
<
0
α
i
=
0
β
i
=
0
0
0
True/false
0
0
0
0
0
0
1
False
ω
+1
−1
0
±1
0
1
True
ω
−1
+1
±1
0
1
0
False
ω
−1
+1
±1
0
1
0
True
ω
+1
−1
0
±1
1
1
True/false
0
0
0
0
0
x
i is the ith bit of the current chromosome, b
i is the ith bit of the best chromosome, Δ
θ
i is the adjusting angle step, and s
(
α
i
,
β
i
) represents the rotating angle direction.In order to balance the search efficiency and search precision, it is necessary to adopt the dynamic adjustment strategy. The dynamic adjustment strategy of the quantum gate based on the expansion coefficient [23] is defined as(9)
ω
=
ω
m
i
n
+
ω
m
a
x
-
ω
m
i
n
1
-
t
i
t
m
a
x
ε
,where ω
m
i
n is the minimum value, ω
m
a
x is the maximum value, t
i is the current genetic generation, t
m
a
x is the largest genetic generation, and ε is the expansion coefficient.


### 3.3. Screw Remaining Life Prediction Model

During the screw accelerated test, accelerometers connected with the online monitoring system are used to monitor ball screw’s vibration. With the method of wavelet packet decomposition, each sensor signal is decomposed to different subbands and original features are extracted. To obtain sensitive feature vectors, Isomap is used to construct the nonlinear feature space mapping from high dimension to low dimension. During the accelerated test, sensitive feature vectors are stored into the database. Meanwhile, motor encoder’s signal reflects the rotating speed of the testing ball screw and torque sensor’s signal reflects the axial force of the testing ball screw. The working condition is composed of these two signals and accumulated in the database (during the accelerated test, the motor speed contains 100 r/min, 300 r/min, and 800 r/min). Obviously, the database consists of sensitive feature vectors and working condition during the whole period of the accelerated test.Samples contained sensitive feature vectors and working condition are randomly chosen from the database. Then the training samples and testing samples can be constructed. In order to get suitable initialization parameters of SVM, training samples are adopted in QGA. In QGA, the training samples are divided into 3 subsamples. Each subsample is used as the validation set while the other two subsamples are used as the training set. Therefore, the prediction accuracy of each subsample can be obtained. Thus, the MSE of these prediction results constitute the fitness function according to (7). By constantly searching and calculating, the best parameters with best fitness can be searched and treated as the optimized SVM parameters. With the optimized SVM parameters, SVM can be initialized. Then the training samples are input into SVM. According to (2)–(4), the SVM is trained by iterative calculation. After that, the testing samples are applied to check the trained SVM. To check the prediction accuracy of the trained SVM, the mean absolute error (MAE) of the prediction result for testing samples is calculated. If the MAE of the prediction result for testing samples is less than 3%, this trained SVM can be used as screw remaining life prediction model. Otherwise the SVM must be retrained with methods of increasing samples or modifying network parameters. Figure 6 shows the flowchart of building screw remaining life prediction model.Figure 6
Flowchart of building screw remaining life prediction model.


## 3.1. Support Vector Machine

In this paper, SVM, which is suitable for studying small sample problems, is taken as the theoretical basis of modeling. As a machine learning algorithm developed from statistical learning theory, SVM maps inseparable learning samples from low-dimensional space into high-dimensional space through a kernel function so as to obtain an optimal hyperplane. Figure2 shows the structure of SVM. If x
i
(
i
=
1,2
,
…
,
N
) are the training samples, y
i are the training goal associated with x
i, and then searching the optimal hyperplane is transformed into solving the convex quadratic equation which can be defined as [20–22](2)
M
i
n
i
m
i
z
e
:
1
2
w
2
Subject  to:
y
i
w
·
x
i
+
b
≥
1
,
i
=
1,2
,
…
,
N
,where w and b are the coefficients of the hyperplane equation.Figure 2
The structure of SVM.According to the solution of convex quadratic programming in the optimization theory, Lagrange function is constructed as(3)
L
w
,
α
,
b
=
1
2
w
2
-
∑
i
=
1
N
α
i
y
i
w
·
x
i
+
b
+
∑
i
=
1
N
α
i
,where α
i is the Lagrange multiplier.To be applied in regression analysis, penalty factorC is introduced in the objective function. Then (2) becomes(4)
M
i
n
i
m
i
z
e
:
∑
i
=
1
N
α
i
-
1
2
∑
i
=
1
N
∑
j
=
1
N
α
i
α
j
y
i
y
j
x
i
·
x
j
K
x
i
,
x
j
Subject  to:
∑
i
=
1
N
α
i
y
j
=
0
,where the range of α
i is [
0
,
C
], i
=
1,2
,
…
,
N.In SVM model, different kernel functions and kernel function parameters have a great impact on the calculation results during the process of mapping the low-dimensional space to high-dimensional space. Meanwhile, the value of penalty factorC also determines the regression effect of SVM. To obtain the best prediction result of SVM, it is significant to choose suitable kernel function parameter and penalty factor C for SVM. Usually, the values of the kernel function parameter and penalty factor C are based on experience. Thus it is difficult to guarantee that the prediction result is the best. Therefore, QGA, a kind of heuristic algorithm, is applied to search the optimal kernel function parameter and penalty factor so as to obtain the best regression prediction effect for SVM model.


## 3.2. Quantum Genetic Algorithm

The main procedure of SVM parameters optimization can be described as follows: according to the variation ranges of kernel function parameter and penalty factor, constructing chromosomal genes with quantum bit coding system, so the evolutionary population which includes several chromosomes can be generated. Then realizing the continuous population evolution by means of cross, variation and quantum rotation gate. Finally, the optimized kernel function parameters and penalty factor can be obtained. Figure3 shows the algorithm flowchart of SVM parameters optimized by QGA.Figure 3
Algorithm flowchart of SVM parameters optimized by QGA.


### 3.2.1. Calculation of Fitness Function

According to the quantum bit coding method, arbitrary chromosome of the populationQ
(
t
0
) can be defined as [23, 24](5)
q
j
t
0
=
α
11
t
0
β
11
t
0
∣
α
12
t
0
β
12
t
0
∣
⋯
⋯
∣
α
1
k
1
t
0
β
1
k
1
t
0
∣
α
21
t
0
β
21
t
0
∣
α
22
t
0
β
22
t
0
∣
⋯
⋯
∣
α
2
k
2
t
0
β
2
k
2
t
0
∣
⋯
⋯
∣
α
m
1
t
0
β
m
1
t
0
∣
α
m
2
t
0
β
m
2
t
0
∣
⋯
⋯
∣
α
m
k
m
t
0
β
m
k
m
t
0
,where q
j
t
0 is thejth chromosome of the t
0 generation and m is the number of the chromosomal genes. The quantum bit number of each chromosomal gene is shown by k
1
,
k
2
,
…
,
k
m.Then population ofN chromosomes can be expressed as(6)
Q
t
0
=
q
1
t
0
,
q
2
t
0
,
…
,
q
N
t
0
.In this paper, SVM kernel function parameter and penalty factor are corresponding to the genes of chromosome q
j
t
0. Since each chromosome should express same initialization state, the initial probability amplitude α
,
β
T is 1
/
2
,
1
/
2
T. As a kind of validation method for assessing the results of a statistical analysis, cross validation is applied to construct the fitness function, so the overfitting problem can be avoided. The key steps of cross validation is to segment a sample of data into several subsets, performing the analysis on one subset which is called the training set, and validating the analysis on the other subset which is called the testing set. The schematic diagram of cross validation is shown in Figure 4.Figure 4
The schematic diagram of cross validation.Then the fitness function can be formed as(7)
f
=
-
1
s
∑
p
=
1
s
M
p
,where s is the number of iterations. M
p represents the mean squared error (MSE) of the prediction results forpth testing set.


### 3.2.2. Quantum Cross and Variation

Quantum cross means exchanging the corresponding genes of any two chromosomes in order to form two similar chromosomes. The purpose of the quantum cross is to exchange information between two possible solutions. The schematic diagram of quantum cross is shown in Figure5.Figure 5
The schematic diagram of quantum cross.By changing one or several genes of any chromosome, a new chromosome can be formed, so as to ensure the diversity of the chromosomes. This process is called quantum variation which is a good way to change the current evolution direction and avoid local optimum.


### 3.2.3. Quantum Rotation Gate

Quantum gate is adopted to be the actuator of the evolution process in QGA. The quantum bit probability amplitude is updated based on quantum rotating angleθ
(8)
α
′
β
′
=
c
o
s
θ
-
s
i
n
θ
s
i
n
θ
c
o
s
θ
α
β
,where α
′
β
′
T is the new quantum bit probability amplitude and α
β
T is the quantum bit probability amplitude before updating.The adjustment strategy [25] of the quantum rotating angle θ is shown in Table 1.Table 1
The adjustment strategy of the quantum rotating angle.
x
i
b
i
f
(
x
)
>
f
(
b
)
Δ
θ
i
s
(
α
i
,
β
i
)
α
i
β
i
>
0
α
i
β
i
<
0
α
i
=
0
β
i
=
0
0
0
True/false
0
0
0
0
0
0
1
False
ω
+1
−1
0
±1
0
1
True
ω
−1
+1
±1
0
1
0
False
ω
−1
+1
±1
0
1
0
True
ω
+1
−1
0
±1
1
1
True/false
0
0
0
0
0
x
i is the ith bit of the current chromosome, b
i is the ith bit of the best chromosome, Δ
θ
i is the adjusting angle step, and s
(
α
i
,
β
i
) represents the rotating angle direction.In order to balance the search efficiency and search precision, it is necessary to adopt the dynamic adjustment strategy. The dynamic adjustment strategy of the quantum gate based on the expansion coefficient [23] is defined as(9)
ω
=
ω
m
i
n
+
ω
m
a
x
-
ω
m
i
n
1
-
t
i
t
m
a
x
ε
,where ω
m
i
n is the minimum value, ω
m
a
x is the maximum value, t
i is the current genetic generation, t
m
a
x is the largest genetic generation, and ε is the expansion coefficient.


## 3.2.1. Calculation of Fitness Function

According to the quantum bit coding method, arbitrary chromosome of the populationQ
(
t
0
) can be defined as [23, 24](5)
q
j
t
0
=
α
11
t
0
β
11
t
0
∣
α
12
t
0
β
12
t
0
∣
⋯
⋯
∣
α
1
k
1
t
0
β
1
k
1
t
0
∣
α
21
t
0
β
21
t
0
∣
α
22
t
0
β
22
t
0
∣
⋯
⋯
∣
α
2
k
2
t
0
β
2
k
2
t
0
∣
⋯
⋯
∣
α
m
1
t
0
β
m
1
t
0
∣
α
m
2
t
0
β
m
2
t
0
∣
⋯
⋯
∣
α
m
k
m
t
0
β
m
k
m
t
0
,where q
j
t
0 is thejth chromosome of the t
0 generation and m is the number of the chromosomal genes. The quantum bit number of each chromosomal gene is shown by k
1
,
k
2
,
…
,
k
m.Then population ofN chromosomes can be expressed as(6)
Q
t
0
=
q
1
t
0
,
q
2
t
0
,
…
,
q
N
t
0
.In this paper, SVM kernel function parameter and penalty factor are corresponding to the genes of chromosome q
j
t
0. Since each chromosome should express same initialization state, the initial probability amplitude α
,
β
T is 1
/
2
,
1
/
2
T. As a kind of validation method for assessing the results of a statistical analysis, cross validation is applied to construct the fitness function, so the overfitting problem can be avoided. The key steps of cross validation is to segment a sample of data into several subsets, performing the analysis on one subset which is called the training set, and validating the analysis on the other subset which is called the testing set. The schematic diagram of cross validation is shown in Figure 4.Figure 4
The schematic diagram of cross validation.Then the fitness function can be formed as(7)
f
=
-
1
s
∑
p
=
1
s
M
p
,where s is the number of iterations. M
p represents the mean squared error (MSE) of the prediction results forpth testing set.


## 3.2.2. Quantum Cross and Variation

Quantum cross means exchanging the corresponding genes of any two chromosomes in order to form two similar chromosomes. The purpose of the quantum cross is to exchange information between two possible solutions. The schematic diagram of quantum cross is shown in Figure5.Figure 5
The schematic diagram of quantum cross.By changing one or several genes of any chromosome, a new chromosome can be formed, so as to ensure the diversity of the chromosomes. This process is called quantum variation which is a good way to change the current evolution direction and avoid local optimum.


## 3.2.3. Quantum Rotation Gate

Quantum gate is adopted to be the actuator of the evolution process in QGA. The quantum bit probability amplitude is updated based on quantum rotating angleθ
(8)
α
′
β
′
=
c
o
s
θ
-
s
i
n
θ
s
i
n
θ
c
o
s
θ
α
β
,where α
′
β
′
T is the new quantum bit probability amplitude and α
β
T is the quantum bit probability amplitude before updating.The adjustment strategy [25] of the quantum rotating angle θ is shown in Table 1.Table 1
The adjustment strategy of the quantum rotating angle.
x
i
b
i
f
(
x
)
>
f
(
b
)
Δ
θ
i
s
(
α
i
,
β
i
)
α
i
β
i
>
0
α
i
β
i
<
0
α
i
=
0
β
i
=
0
0
0
True/false
0
0
0
0
0
0
1
False
ω
+1
−1
0
±1
0
1
True
ω
−1
+1
±1
0
1
0
False
ω
−1
+1
±1
0
1
0
True
ω
+1
−1
0
±1
1
1
True/false
0
0
0
0
0
x
i is the ith bit of the current chromosome, b
i is the ith bit of the best chromosome, Δ
θ
i is the adjusting angle step, and s
(
α
i
,
β
i
) represents the rotating angle direction.In order to balance the search efficiency and search precision, it is necessary to adopt the dynamic adjustment strategy. The dynamic adjustment strategy of the quantum gate based on the expansion coefficient [23] is defined as(9)
ω
=
ω
m
i
n
+
ω
m
a
x
-
ω
m
i
n
1
-
t
i
t
m
a
x
ε
,where ω
m
i
n is the minimum value, ω
m
a
x is the maximum value, t
i is the current genetic generation, t
m
a
x is the largest genetic generation, and ε is the expansion coefficient.


## 3.3. Screw Remaining Life Prediction Model

During the screw accelerated test, accelerometers connected with the online monitoring system are used to monitor ball screw’s vibration. With the method of wavelet packet decomposition, each sensor signal is decomposed to different subbands and original features are extracted. To obtain sensitive feature vectors, Isomap is used to construct the nonlinear feature space mapping from high dimension to low dimension. During the accelerated test, sensitive feature vectors are stored into the database. Meanwhile, motor encoder’s signal reflects the rotating speed of the testing ball screw and torque sensor’s signal reflects the axial force of the testing ball screw. The working condition is composed of these two signals and accumulated in the database (during the accelerated test, the motor speed contains 100 r/min, 300 r/min, and 800 r/min). Obviously, the database consists of sensitive feature vectors and working condition during the whole period of the accelerated test.Samples contained sensitive feature vectors and working condition are randomly chosen from the database. Then the training samples and testing samples can be constructed. In order to get suitable initialization parameters of SVM, training samples are adopted in QGA. In QGA, the training samples are divided into 3 subsamples. Each subsample is used as the validation set while the other two subsamples are used as the training set. Therefore, the prediction accuracy of each subsample can be obtained. Thus, the MSE of these prediction results constitute the fitness function according to (7). By constantly searching and calculating, the best parameters with best fitness can be searched and treated as the optimized SVM parameters. With the optimized SVM parameters, SVM can be initialized. Then the training samples are input into SVM. According to (2)–(4), the SVM is trained by iterative calculation. After that, the testing samples are applied to check the trained SVM. To check the prediction accuracy of the trained SVM, the mean absolute error (MAE) of the prediction result for testing samples is calculated. If the MAE of the prediction result for testing samples is less than 3%, this trained SVM can be used as screw remaining life prediction model. Otherwise the SVM must be retrained with methods of increasing samples or modifying network parameters. Figure 6 shows the flowchart of building screw remaining life prediction model.Figure 6
Flowchart of building screw remaining life prediction model.


## 4. The Experiment Results

### 4.1. Signal Processing and Sensitive Feature Vectors

In order to get a fast calculation speed, Harr wavelet function is chosen as the wavelet packet decomposition function. By wavelet packet decomposition method, each accelerometer signal is decomposed into 5 levels, so 32 features can be obtained in each direction. In accelerated test, the time interval of data acquisition is 20 hours, which we call as a degradation cycle. Therefore, the whole life cycle is divided into 958 degradation cycles considering the rated life of the testing ball screw is approximately 19167 hours. Figures7(a), 7(b), and 7(c) show the variations of original features in three directions. It is clear that the number of original features is 96 in total. If so many original features are directly input into the prediction model, it is easy to cause the dimensionality disaster of the prediction model, which may lead to that the prediction model cannot be trained. In order to prove the validity of the sensitive features selected by Isomap, principal component analysis (PCA) is also introduced to compare with Isomap. PCA is a kind of dimensionality reduction algorithm which is widely adopted in machinery and equipment fault diagnosis. Figure 8(a) shows the sensitive features obtained by Isomap while Figure 8(b) shows the sensitive features obtained by PCA. We can see that the number of the original features are both reduced to 19 through Isomap or PCA, thus these two methods are both succeeded in reducing the dimensionality of the original features. To compare these two kinds of sensitive features, correlation coefficient between sensitive features and the degradation cycle is defined as(10)
r
x
y
=
∑
k
=
1
n
x
k
-
x
-
y
k
-
y
-
∑
k
=
1
n
x
k
-
x
-
2
∑
k
=
1
n
y
k
-
y
-
2
,where x
k is a sensitive feature vector and y
k is the degradation cycle.Figure 7
Normalized values of original features: (a) features ofX direction, (b) features of Y direction, and (c) features of Z direction.
(a)
(b)
(c)Figure 8
Normalized values of sensitive features: (a) features after Isomap and (b) features after PCA.
(a)
(b)Admittedly, correlation coefficient between sensitive feature and the degradation cycle can be used to evaluate feature and the more the better. The comparison between Isomap and PCA is shown in Table2. It is clear that the sensitive features obtained by Isomap contain better tendency with degradation cycle than the sensitive features obtained by PCA. Thus the effectiveness of the Isomap algorithm is demonstrated.Table 2
Accumulation of the correlation coefficients between sensitive features and degradation cycle.
Accumulation of the correlation coefficients
The first
The first five
The first ten
The first fifteen
The first nineteen
Isomap
0.7314
2.1017
2.9849
3.7198
3.9438
PCA
0.7029
1.8792
2.6720
3.2285
3.6775


### 4.2. Screw Remaining Life Prediction

For ball screw, the rate lifeL
h can be defined as(11)
L
h
=
1
60
×
n
m
C
a
F
m
3
×
1
0
6
,where n
m is the rotational speed of ball screw,
C
a is the rated dynamic load, and F
m is the axial load.Thus(12)
L
e
=
L
h
-
L
w
L
h
,where L
e is the normalized remaining life and L
w is the servicing life.In order to validate the effect of the proposed method, back propagation neural network (BP, a kind of feedforward neural network) and Elman neural network (a kind of feedback neural network) are also adopted to train and test with the same samples. For BP neural network, the number of iterations is 100 and the learning rate is 0.1. Meanwhile, the number of iterations is 1000 for Elman neural network. In order to reduce the impact of stochastic factors on BP and Elman neural networks, each model would be built 20 times. Then the average performance is taken as the final performance of the neural network.Moreover, SVM with random parameters and experiential parameters are also applied to compare. 200 samples chosen from 100 degradation cycles are randomly selected from the database; thus training samples include 100 samples and testing samples include 100 samples. For QGA, the population scale and the number of iterations should be based on comprehensive consideration of computational efficiency and sample size. In this paper, initial parameters of QGA are chosen as follows: population scale is 60, the number of iterations is 50, cross probability is 0.3, variation probability is 0.1, and expansion coefficient is 2. The searching interval of SVM parameters is[
-
8,8
], which means that searching intervals of log
2
⁡
c and log
2
⁡
g are both belong to [
-
8,8
]. Prediction models established by BP, Elman, SVM with random parameters, SVM with experiential parameters and SVM with parameters based on QGA are trained with 100 samples chosen from 100 degradation cycles. Figures 9, 10, 11, 12, and 13 show the testing results of five models with another 100 testing samples. Among them, Figures 9(a), 10(a), 11(a), 12(a), and 13(a) show the prediction curves while Figures 9(b), 10(b), 11(b), 12(b), and 13(b) show the prediction errors. It can be seen that models established by SVM with experiential parameters and parameters based on QGA got better prediction accuracy than other models. To check out which parameters are the best choice, we calculate the MSE and MAE of the testing samples. The comparisons of five models are shown in Table 3.Table 3
Comparison of five methods for predicting.
log
2
⁡
c
log
2
⁡
g
MSE of cross validation
MSE of testing samples
MAE of testing samples
BP
~
~
~
4.476 × 10−4
0.0175
Elman
~
~
~
4.638 × 10−4
0.0176
SVM with random parameters
−5
0
0.0428
5.633 × 10−4
0.0201
SVM with experiential parameters
0
−3
0.0174
3.426 × 10−4
0.0158
SVM with parameters based on QGA
0.0078
−2.479
0.0171
2.816 × 10−4
0.0144Figure 9
Prediction results of BP: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 10
Prediction results of Elman: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 11
Prediction results of the SVM with random parameters: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 12
Prediction results of the SVM with experiential parameters: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 13
Prediction results of the SVM with parameters based on QGA: (a) prediction curve and (b) prediction error.
(a)
(b)From Table3, we can know that model established by SVM with parameters based on QGA achieves best prediction effect since both MSE and MAE of the testing samples are minimum. MSE of cross validation, which is the opposite number of the fitness function in QGA, reflects a very good consistency with the MSE and MAE of the testing samples though it is a kind of calculation result from the training samples. Therefore, it is clear that cross validation method is a good way to avoid overtraining and can be used to construct fitness function in QGA.To better observe MSE distribution of cross validation, method of discrete searching is applied; thuslog
2
⁡
c and log
2
⁡
g are both discretized within a range of [
-
8,8
]. By calculating the MSE on each discrete point, a MSE distribution with a range of [
-
8,8
] can be obtained. Figure 14 shows the MSE distribution. It can be seen that the MSE distribution plane is rugged, while experiential parameters fall in the area of less error, but experiential parameters are usually not the optimal parameters of the minimum error. Combined with the results presented in Table 3, it is clear that QGA can effectively search the best parameters of SVM so the screw remaining life prediction model achieves a good prediction accuracy.Figure 14
MSE distribution of cross validation.To test the stability of the screw remaining life prediction model, data of 101st degradation cycle to 150th degradation cycle are adopted. Figure15 shows the prediction results of the proposed method. It is clear that the model maintains a good prediction accuracy for the next 50 degradation cycles. It is shown that the proposed method simultaneously exhibits good stability and precision in processing of the performance degradation data of ball screw. As is shown in Table 4, we can see that the proposed model achieves best prediction effect. The MAE of predicting samples is only 1.11%, which is the smallest among five methods.Table 4
Prediction results of 101st degradation cycle to 150th degradation cycle.
BP
Elman
SVM with random parameters
SVM with experiential parameters
SVM with parameters based on QGA
MSE of predicting samples
3.200 × 10−3
3.000 × 10−3
2.412 × 10−4
2.196 × 10−4
2.102 × 10−4
MAE of predicting samples
0.0522
0.0492
0.0127
0.0113
0.0111Figure 15
Prediction results of 101st degradation cycle to 150th degradation cycle: (a) prediction curve and (b) prediction error.
(a)
(b)


## 4.1. Signal Processing and Sensitive Feature Vectors

In order to get a fast calculation speed, Harr wavelet function is chosen as the wavelet packet decomposition function. By wavelet packet decomposition method, each accelerometer signal is decomposed into 5 levels, so 32 features can be obtained in each direction. In accelerated test, the time interval of data acquisition is 20 hours, which we call as a degradation cycle. Therefore, the whole life cycle is divided into 958 degradation cycles considering the rated life of the testing ball screw is approximately 19167 hours. Figures7(a), 7(b), and 7(c) show the variations of original features in three directions. It is clear that the number of original features is 96 in total. If so many original features are directly input into the prediction model, it is easy to cause the dimensionality disaster of the prediction model, which may lead to that the prediction model cannot be trained. In order to prove the validity of the sensitive features selected by Isomap, principal component analysis (PCA) is also introduced to compare with Isomap. PCA is a kind of dimensionality reduction algorithm which is widely adopted in machinery and equipment fault diagnosis. Figure 8(a) shows the sensitive features obtained by Isomap while Figure 8(b) shows the sensitive features obtained by PCA. We can see that the number of the original features are both reduced to 19 through Isomap or PCA, thus these two methods are both succeeded in reducing the dimensionality of the original features. To compare these two kinds of sensitive features, correlation coefficient between sensitive features and the degradation cycle is defined as(10)
r
x
y
=
∑
k
=
1
n
x
k
-
x
-
y
k
-
y
-
∑
k
=
1
n
x
k
-
x
-
2
∑
k
=
1
n
y
k
-
y
-
2
,where x
k is a sensitive feature vector and y
k is the degradation cycle.Figure 7
Normalized values of original features: (a) features ofX direction, (b) features of Y direction, and (c) features of Z direction.
(a)
(b)
(c)Figure 8
Normalized values of sensitive features: (a) features after Isomap and (b) features after PCA.
(a)
(b)Admittedly, correlation coefficient between sensitive feature and the degradation cycle can be used to evaluate feature and the more the better. The comparison between Isomap and PCA is shown in Table2. It is clear that the sensitive features obtained by Isomap contain better tendency with degradation cycle than the sensitive features obtained by PCA. Thus the effectiveness of the Isomap algorithm is demonstrated.Table 2
Accumulation of the correlation coefficients between sensitive features and degradation cycle.
Accumulation of the correlation coefficients
The first
The first five
The first ten
The first fifteen
The first nineteen
Isomap
0.7314
2.1017
2.9849
3.7198
3.9438
PCA
0.7029
1.8792
2.6720
3.2285
3.6775


## 4.2. Screw Remaining Life Prediction

For ball screw, the rate lifeL
h can be defined as(11)
L
h
=
1
60
×
n
m
C
a
F
m
3
×
1
0
6
,where n
m is the rotational speed of ball screw,
C
a is the rated dynamic load, and F
m is the axial load.Thus(12)
L
e
=
L
h
-
L
w
L
h
,where L
e is the normalized remaining life and L
w is the servicing life.In order to validate the effect of the proposed method, back propagation neural network (BP, a kind of feedforward neural network) and Elman neural network (a kind of feedback neural network) are also adopted to train and test with the same samples. For BP neural network, the number of iterations is 100 and the learning rate is 0.1. Meanwhile, the number of iterations is 1000 for Elman neural network. In order to reduce the impact of stochastic factors on BP and Elman neural networks, each model would be built 20 times. Then the average performance is taken as the final performance of the neural network.Moreover, SVM with random parameters and experiential parameters are also applied to compare. 200 samples chosen from 100 degradation cycles are randomly selected from the database; thus training samples include 100 samples and testing samples include 100 samples. For QGA, the population scale and the number of iterations should be based on comprehensive consideration of computational efficiency and sample size. In this paper, initial parameters of QGA are chosen as follows: population scale is 60, the number of iterations is 50, cross probability is 0.3, variation probability is 0.1, and expansion coefficient is 2. The searching interval of SVM parameters is[
-
8,8
], which means that searching intervals of log
2
⁡
c and log
2
⁡
g are both belong to [
-
8,8
]. Prediction models established by BP, Elman, SVM with random parameters, SVM with experiential parameters and SVM with parameters based on QGA are trained with 100 samples chosen from 100 degradation cycles. Figures 9, 10, 11, 12, and 13 show the testing results of five models with another 100 testing samples. Among them, Figures 9(a), 10(a), 11(a), 12(a), and 13(a) show the prediction curves while Figures 9(b), 10(b), 11(b), 12(b), and 13(b) show the prediction errors. It can be seen that models established by SVM with experiential parameters and parameters based on QGA got better prediction accuracy than other models. To check out which parameters are the best choice, we calculate the MSE and MAE of the testing samples. The comparisons of five models are shown in Table 3.Table 3
Comparison of five methods for predicting.
log
2
⁡
c
log
2
⁡
g
MSE of cross validation
MSE of testing samples
MAE of testing samples
BP
~
~
~
4.476 × 10−4
0.0175
Elman
~
~
~
4.638 × 10−4
0.0176
SVM with random parameters
−5
0
0.0428
5.633 × 10−4
0.0201
SVM with experiential parameters
0
−3
0.0174
3.426 × 10−4
0.0158
SVM with parameters based on QGA
0.0078
−2.479
0.0171
2.816 × 10−4
0.0144Figure 9
Prediction results of BP: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 10
Prediction results of Elman: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 11
Prediction results of the SVM with random parameters: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 12
Prediction results of the SVM with experiential parameters: (a) prediction curve and (b) prediction error.
(a)
(b)Figure 13
Prediction results of the SVM with parameters based on QGA: (a) prediction curve and (b) prediction error.
(a)
(b)From Table3, we can know that model established by SVM with parameters based on QGA achieves best prediction effect since both MSE and MAE of the testing samples are minimum. MSE of cross validation, which is the opposite number of the fitness function in QGA, reflects a very good consistency with the MSE and MAE of the testing samples though it is a kind of calculation result from the training samples. Therefore, it is clear that cross validation method is a good way to avoid overtraining and can be used to construct fitness function in QGA.To better observe MSE distribution of cross validation, method of discrete searching is applied; thuslog
2
⁡
c and log
2
⁡
g are both discretized within a range of [
-
8,8
]. By calculating the MSE on each discrete point, a MSE distribution with a range of [
-
8,8
] can be obtained. Figure 14 shows the MSE distribution. It can be seen that the MSE distribution plane is rugged, while experiential parameters fall in the area of less error, but experiential parameters are usually not the optimal parameters of the minimum error. Combined with the results presented in Table 3, it is clear that QGA can effectively search the best parameters of SVM so the screw remaining life prediction model achieves a good prediction accuracy.Figure 14
MSE distribution of cross validation.To test the stability of the screw remaining life prediction model, data of 101st degradation cycle to 150th degradation cycle are adopted. Figure15 shows the prediction results of the proposed method. It is clear that the model maintains a good prediction accuracy for the next 50 degradation cycles. It is shown that the proposed method simultaneously exhibits good stability and precision in processing of the performance degradation data of ball screw. As is shown in Table 4, we can see that the proposed model achieves best prediction effect. The MAE of predicting samples is only 1.11%, which is the smallest among five methods.Table 4
Prediction results of 101st degradation cycle to 150th degradation cycle.
BP
Elman
SVM with random parameters
SVM with experiential parameters
SVM with parameters based on QGA
MSE of predicting samples
3.200 × 10−3
3.000 × 10−3
2.412 × 10−4
2.196 × 10−4
2.102 × 10−4
MAE of predicting samples
0.0522
0.0492
0.0127
0.0113
0.0111Figure 15
Prediction results of 101st degradation cycle to 150th degradation cycle: (a) prediction curve and (b) prediction error.
(a)
(b)


## 5. Conclusion

In this paper, screw remaining life prediction method based on QGA and SVM is proposed. The experiment results show that our proposed method achieves best prediction accuracy compared to another two models. Some conclusions can be got as follows:(1)
Acceleration signals of screw nut can be acquired to monitor the performance degradation of ball screw.
(2)
Combined wavelet packet decomposition and Isomap method, the sensitive feature vectors can be extracted.
(3)
To improve the prediction accuracy of SVM model, an optimization algorithm based on QGA and cross validation is presented in this paper. Compared with models adopting random parameters and experiential parameters, model with the optimal parameters achieves the best prediction accuracy. Therefore, proposed approach is suitable for screw remaining life prediction.



---

*Source: 9581379-2017-02-15.xml*