Ruthger Righart
Email: rrighart@googlemail.com
Blog website: https://rrighart.github.io
Hire-me remotely, for short / long projects: https://rranalytics.weebly.com/
1. Introduction
Can one predict when an engine or device breaks down?
This seems a pure engineering question. But nowadays it is also a data science question. More concretely, it is an important question everywhere engines and devices use data to guide maintenance, such as aircraft engines1, windturbine engines2, and rotating machinery3. With regards to human safety and logistic planning, it is not a good idea to just wait until an engine breaks down. It is essential to proactively plan maintenance to avoid costly downtime.
Proactive planning is enabled by multiple sensor technologies that provide powerful data about the state of machines and devices. Sensors include data such as temperature, fan and core speed, static pressure etc. Can we use these data to predict within certain margins how long an aircraft engine will be functioning without failure? And if so, how to do it?
This is the question the concept of remaining useful life (RUL) attempts to answer. It aims to estimate the remaining time an item, component, or system is able to function in accordance with its intended purpose before warranting replacement. The present blog shows how to use deep learning in Python Keras to predict the RUL. It is meant to provide an example case study, not an exhaustive and ultimate solution.
There is a lack of real data to answer this question. However, data simulations have been made and provide a unique resource. One such a fascinating simulation is provided by the C-MAPSS data1. It provides train data that show sensor-based time-series until the timepoint the engine breaks down. In contrast, the test data constitute of sensor-based time-series a "random" time before the endpoint. The key aim is to estimate the RUL of the testset, that is, how much time is left after the last recorded time-point.
2. Modules
Python 2 was used for the current blog. Python 3 can be used with a few adaptations to the code. The following modules will be used:
import os
import pandas as pd
import numpy as np
np.random.seed(1337)
import requests, zipfile, StringIO
from IPython.display import Image
import matplotlib as mpl
import matplotlib.pyplot as plt
from pandas import read_csv
from matplotlib import pyplot as plt
from sklearn.metrics import mean_squared_error
from keras.models import Sequential
from keras.layers.core import Dense, Activation
from keras.wrappers.scikit_learn import KerasRegressor
from sklearn.model_selection import cross_val_score
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import MinMaxScaler
The following settings will be used to avoid exponential values in output or tables and to display 50 rows maximum:
pd.set_option('display.float_format', lambda x: '%.3f' % x)
pd.options.display.max_rows=50
3. Sketch of the question
A simplified lineplot illustrates best the question at hand. Given a fictive temperature sensor, for which we have 8 cycles, are we able to predict the remaining cycles?
The kind of questions that can be addressed for such a time-serie:
- Can we forecast the temperature of the future time-points?
- Can we predict how long the time-serie continues until a pre-defined event happens?
A=[60,63,67,74,77,81,82,87,92]
B=[92,99,104,110,116,125]
C = np.append(np.repeat(np.nan, len(A)-1), B)
plt.figure(figsize = (16, 8))
plt.plot(A, color='red', linewidth=3)
plt.plot(C, 'r:', linewidth=3)
plt.axvline(x=len(A)-1, color='grey', linestyle=':', linewidth=4)
plt.axvline(x=len(C)-1, color='grey', linestyle=':', linewidth=4)
plt.title('Example temperature sensor', fontsize=16)
plt.xlabel('# Cycles', fontsize=16)
plt.ylabel('Degrees', fontsize=16)
plt.xticks(fontsize=16)
plt.yticks(fontsize=16)
plt.show()
4. Loading the data
The data concerns the C-MAPSS set and can be found at: https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/ where it is under "Turbofan Engine Degradation Simulation Data Set" providing a ZIP file. The one that we are using for the current purpose is FD001, but the reader is encouraged to use the others as well (more data usually means better prediction power). It can be downloaded either manually or using the script below that downloads the content of the ZIP file from the requested URL in your homedirectory4 .
r = requests.get('https://ti.arc.nasa.gov/c/6/', stream=True)
z = zipfile.ZipFile(StringIO.StringIO(r.content))
z.extractall()
Now that the files are in your home directory, you can read them using Pandas read_csv. We will load a train- and testset, as well as a RUL set. RUL contains the true values for remaining useful life to which we are going to compare our predicted values in the testset.
train = pd.read_csv('train_FD001.csv', parse_dates=False, delimiter=" ", decimal=".", header=None)
test = pd.read_csv('test_FD001.csv', parse_dates=False, delimiter=" ", decimal=".", header=None)
RUL = pd.read_csv('RUL_FD001.csv', parse_dates=False, delimiter=" ", decimal=".", header=None)
5. Missing values
First, we need to drop the two last columns since they actually only have missing values. This is probably due to trailing tab characters in the csv file. At any time, it is always better to verify this as the owners of the data may have edited this at any moment. The following code prints the proportion of missing values for each column in the train- and testset:
tableNA = pd.concat([train.isnull().sum(), test.isnull().sum()], axis=1)
tableNA.columns = ['train', 'test']
tableNA
We will first drop the columns that consisted of missing values:
train.drop(train.columns[[-1,-2]], axis=1, inplace=True)
test.drop(test.columns[[-1,-2]], axis=1, inplace=True)
And then we will give names to all columns:
cols = ['unit', 'cycles', 'op_setting1', 'op_setting2', 'op_setting3', 's1', 's2', 's3', 's4', 's5', 's6', 's7', 's8', 's9', 's10', 's11', 's12', 's13', 's14', 's15', 's16', 's17', 's18', 's19', 's20', 's21']
train.columns = cols
test.columns = cols
And so the train- and testset look as follows:
train.head()
test.head()
And the RUL data look as follows, in short meaning that the remaining useful life for the first unit was 112 cycles, the second unit 98 cycles, etc.
RUL.head()
6. Outliers and flat lines
To know if there are outliers (extreme values) in the data, we could use descriptive statistics, train.describe().transpose(), and see if the min. and max. values are far away from the central tendency. As we can see below, this is not the case for any of the sensors.
However, if we look carefully we can see something else that is quite remarkable: there are several sensors where the min. and max. values are identical, and where the standard deviation (std) is zero. In time-series, this is called a flat line, which means there is no activity, possibly caused by sensor malfunctioning.
train.describe().transpose()
The sensors s1, s5, s10, s16, s18, and s19 as well as op_setting 3, will for this reason be removed from further analyses:
train.drop(['s1', 's5', 's10', 's16', 's18', 's19', 'op_setting3'], axis=1, inplace=True)
test.drop(['s1', 's5', 's10', 's16', 's18', 's19', 'op_setting3'], axis=1, inplace=True)
The distribution of the columns looks as follows:
train.hist(bins=50, figsize=(18,16))
plt.show()
7. Exploratory analyses of the max. number of cycles per unit
Exploratory data analyses provide insight into the aircraft engines in action. For example, it would be good to have an idea of the maximum lifetime of the 100 different units. The barplots below show that there is a large variation across units regarding max. number of cycles, and that, as expected, the number of cycles is shorter for testset than trainset.
cyclestrain = train.groupby('unit', as_index=False)['cycles'].max()
cyclestest = test.groupby('unit', as_index=False)['cycles'].max()
fig = plt.figure(figsize = (16,12))
fig.add_subplot(1,2,1)
bar_labels = list(cyclestrain['unit'])
bars = plt.bar(list(cyclestrain['unit']), cyclestrain['cycles'], color='red')
plt.ylim([0, 400])
plt.xlabel('Units', fontsize=16)
plt.ylabel('Max. Cycles', fontsize=16)
plt.title('Max. Cycles per unit in trainset', fontsize=16)
plt.xticks(np.arange(min(bar_labels)-1, max(bar_labels)-1, 5.0), fontsize=12)
plt.yticks(fontsize=12)
fig.add_subplot(1,2,2)
bars = plt.bar(list(cyclestest['unit']), cyclestest['cycles'], color='grey')
plt.ylim([0, 400])
plt.xlabel('Units', fontsize=16)
plt.ylabel('Max. Cycles', fontsize=16)
plt.title('Max. Cycles per unit in testset', fontsize=16)
plt.xticks(np.arange(min(bar_labels)-1, max(bar_labels)-1, 5.0), fontsize=12)
plt.yticks(fontsize=12)
plt.show()
In fact, in machine learning it is considered good behavior to put apart the testset. So from now we do not touch or look at it anymore.
8. Visualization of several sensors, for a particular unit
The following visualizes different sensors of a particular unit. It gives a good impression that all sensors have a different range of values (as they measure different entities as speed, temperature), and that they do not all show an identical pattern or trend. In fact, some do increase in amplitude, while others decrease in amplitude over time.
values = train[train.unit==1].values
groups = [5, 6, 7, 8, 9, 10, 11,12,13]
i = 1
plt.figure(figsize=(10,20))
for group in groups:
plt.subplot(len(groups), 1, i)
plt.plot(values[:, group])
plt.title(train.columns[group], y=0.5, loc='right')
i += 1
plt.show()
We can also show a single sensor for different engine units. This nicely illustrates across units that amplitudes decrease over time and seem to go to a certain minimum threshold of about 551-552. By the way: the data seem rather noisy and filtering may help at some point (not treated in this blog).
plt.figure(figsize = (8, 8))
plt.plot(train[train.unit==1].cycles, train[train.unit==1].s7)
plt.plot(train[train.unit==2].cycles, train[train.unit==2].s7)
plt.plot(train[train.unit==3].cycles, train[train.unit==3].s7)
plt.plot(train[train.unit==4].cycles, train[train.unit==4].s7)
plt.plot(train[train.unit==5].cycles, train[train.unit==5].s7)
plt.plot(train[train.unit==6].cycles, train[train.unit==6].s7)
plt.plot(train[train.unit==7].cycles, train[train.unit==7].s7)
plt.plot(train[train.unit==8].cycles, train[train.unit==8].s7)
plt.plot(train[train.unit==9].cycles, train[train.unit==9].s7)
plt.plot(train[train.unit==10].cycles, train[train.unit==10].s7)
plt.plot(train[train.unit==11].cycles, train[train.unit==11].s7)
plt.xlabel('# Cycles')
plt.ylabel('Sensor measurements')
plt.show()
So would it be possible that different units have a similar minimum and maximum value for single sensors? This would make sense if sensors started at a low amplitude (for ex. temperature) and went up to a high amplitude over time (or the other way around for other metrics). We could test this for ten units:
minb = train.groupby('unit', as_index=False).min().head(10)
maxb = train.groupby('unit', as_index=False).max().head(10)
mmtable = minb.append(maxb, ignore_index=True)
The following plot suggests that sensors follow a similar kind of pattern for different units.
plt.figure(figsize = (12,12))
col = np.concatenate((np.repeat('red', 10), np.repeat('blue', 10)), axis=0)
bar_labels = list(mmtable['unit'])
x_pos = list(range(len(bar_labels)))
bars = plt.bar(x_pos, mmtable['s2'], color=col)
plt.ylim([640, 645])
plt.xlabel('Units', fontsize=14)
plt.ylabel('Measure s2', fontsize=14)
plt.xticks(x_pos, bar_labels, fontsize=14)
plt.yticks(fontsize=14)
plt.show()
This were some exploratory analyses and visualizations; certainly not exhaustive but a starting point to get insight into the characteristics of the data.
9. Establishing remaining life in cycles
It is now about time to determine the remaining useful life (RUL) for the trainset, for each row. First, we determine in the trainset for each row the max. cycles for the particular unit. We use the groupby command to obtain for every unit the max. number of cycles, and in turn use pd.merge to bring these values into the original train set:
train = pd.merge(train, train.groupby('unit', as_index=False)['cycles'].max(), how='left', on='unit')
train.rename(columns={"cycles_x": "cycles", "cycles_y": "maxcycles"}, inplace=True)
We then determine the time to failure (TTF) for every row, which is the number of cycles subtracted from the maximum number of cycles in a particular unit.
train['TTF'] = train['maxcycles'] - train['cycles']
10. Scaling
Another preparatory step that is important is scaling. We are going to use the MinMaxScaler in Python:
scaler = MinMaxScaler()
Before scaling, let us inspect the original descriptive statistics. This shows that there are huge differences between multiple sensors in minimum and maximum values, as expected since the sensors measure different entities (such as temperature, speed):
train.describe().transpose()
We first make a copy of the data, sothat we have a dataset for unscaled and scaled data:
ntrain = train.copy()
And then we select the data that we would like to scale:
ntrain.iloc[:,2:19] = scaler.fit_transform(ntrain.iloc[:,2:19])
To inspect the scaling, it would be important to see the minimum and maximum value for each column.
ntrain.describe().transpose()
We are going to scale the testdata using the scaler settings of the traindata.
ntest = test.copy()
It concerns the following columns:
pd.DataFrame(ntest.columns).transpose()
ntest.iloc[:,2:19] = scaler.transform(ntest.iloc[:,2:19])
ntest.describe().transpose()
11. Visualize the scaled data
It is always a good idea to visualize the scaled data. This to ensure that the data look similar after scaling (except from the numbers on the Y-axis, of course).
fig = plt.figure(figsize = (8, 8))
fig.add_subplot(1,2,1)
plt.plot(train[train.unit==1].s2)
plt.plot(test[test.unit==1].s2)
plt.legend(['Train','Test'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylabel('Original unit')
fig.add_subplot(1,2,2)
plt.plot(ntrain[ntrain.unit==1].s2)
plt.plot(ntest[ntest.unit==1].s2)
plt.legend(['Scaled Train','Scaled Test'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylabel('Scaled unit')
plt.show()
12. Fraction time to failure
Time to failure (TTF) for each unit has a different length, and it would be good to express this in a fraction of remaining number of cycles. This starts for a particular unit at 1.00, and goes to 0.00, the point where the engine fails (TTFx). It is in fact similar to scaling, but here it is applied at the unit level.
In Python, we can express this in a function:
def fractionTTF(dat,q):
return(dat.TTF[q]-dat.TTF.min()) / float(dat.TTF.max()-dat.TTF.min())
fTTFz = []
fTTF = []
for i in range(train['unit'].min(),train['unit'].max()+1):
dat=train[train.unit==i]
dat = dat.reset_index(drop=True)
for q in range(len(dat)):
fTTFz = fractionTTF(dat, q)
fTTF.append(fTTFz)
ntrain['fTTF'] = fTTF
The following plot shows on the left TTF in cycles, for the first 4 units. On the right, the fraction TTF that starts at 1.00 and ends at 0.00. However, the number of cycles that it takes to failure remains identical. Some units have a longer duration than others, as can be clearly seen in the X-axes of the zigzag figures5
mx = cyclestrain.iloc[0:4,1].sum()
fig = plt.figure(figsize = (8, 8))
fig.add_subplot(1,2,1)
plt.plot(ntrain.TTF[0:mx])
plt.legend(['Time to failure (in cycles)'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylabel('Original unit')
fig.add_subplot(1,2,2)
plt.plot(ntrain.fTTF[0:mx])
plt.legend(['Time to failure (fraction)'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylabel('Scaled unit')
plt.show()
ntrain['fTTF'].describe()
Finally, we have the following columns that we select from for training the model.
pd.DataFrame(ntrain.columns).transpose()
13. Neural network in Keras
We are now ready to train the data and predict the RUL. For this we are using Keras. Importantly, for the target variable Y_train, we take the fraction of the time to failure TTFx. The features are the scaled values. Note that the data are in Numpy.ndarray format, which can be checked by type(X_train)
X_train = ntrain.values[:,1:19]
Y_train = ntrain.values[:, 21]
X_test = ntest.values[:,1:19]
About the neural network in Keras: A ReLU activation function is used. We use the Adam optimizer. The model is not optimized, and so probably still better results can be obtained by tuning the parameters. The neural network that we use here has 18 input neurons, 6 intermediate neurons, 1 output neuron (1 output since it is a regression problem, estimating the remaining useful lifetime).
url = 'https://raw.githubusercontent.com/RRighart/Gatu/master/nn-18-6-1.png'
Image(url)
Everything is ready to train the model:
model = Sequential()
model.add(Dense(6, input_dim=18, kernel_initializer='normal', activation='relu'))
model.add(Dense(1, kernel_initializer='normal'))
model.compile(loss='mean_squared_error', optimizer='adam')
model.fit(X_train, Y_train, nb_epoch=20)
14. Predict test score, fraction RUL from 1.00-0.00
We can now predict the fraction RUL values for the testset:
score = model.predict(X_test)
Remind that the values Y_train were before transformed to a fraction of remaining useful life, from 1.0 till 0.00, in order to cancel out the possible effect of total cycle duration. So score should also contain values in this range:
score[0:10]
The values are roughly in the range 0-1:
print(score.min(), score.max())
Finally, we would like to re-convert the fraction of remaining life (0.00-1.00) to remaining useful life as expressed in number of cycles.
For this, we would first need a column with the maximum number of cycles per unit that were in the testset, which can be obtained by the groupby and merge functions (we did this before for the train set as well):
test = pd.merge(test, test.groupby('unit', as_index=False)['cycles'].max(), how='left', on='unit')
test.rename(columns={"cycles_x": "cycles", "cycles_y": "maxcycles"}, inplace=True)
test['score'] = score
test.head()
Remind that "test" contains only the unscaled values, whereas for modeling and predicting "score" the scaled features were used.
Second, knowing the predicted remaining life (fraction), we need to estimate the predicted total number of cycles per unit in the testset. This can be done with the following function:
def totcycles(data):
return(data['cycles'] / (1-data['score']))
test['maxpredcycles'] = totcycles(test)
Last, we subtract the maximum cycles per unit from the predicted total number of cycles in the testset to obtain the RUL, remaining useful lifetime:
def RULfunction(data):
return(data['maxpredcycles'] - data['maxcycles'])
test['RUL'] = RULfunction(test)
From this column RUL in the testdata we can reconstruct the remaining useful life at the point the maximum cycle is reached, in each unit:
test['RUL'].head()
15. Predict RUL in cycles
The following will compute the RUL per unit (based on the max. cycles) from the RUL column that contains predicted values for each row.
t = test.columns == 'RUL'
ind = [i for i, x in enumerate(t) if x]
predictedRUL = []
for i in range(test.unit.min(), test.unit.max()+1):
npredictedRUL=test[test.unit==i].iloc[test[test.unit==i].cycles.max()-1,ind]
predictedRUL.append(npredictedRUL)
predictedRUL[0:10]
len(predictedRUL)
We had the so-called "zigzag" figure for the trainset. Now we can re-construct and visualize this for the testset, for the predicted and true values. Let us first create a list of values that discounts for every RUL per unit:
xtrueRUL = list(RUL.loc[:,0])
otrueRUL = []
for i in range(0,len(xtrueRUL)):
otrueRUL = np.concatenate((otrueRUL, list(reversed(np.arange(xtrueRUL[i])))))
xpredictedRUL = list(round(x) for x in predictedRUL)
opredictedRUL = []
for i in range(0,len(xpredictedRUL)):
opredictedRUL = np.concatenate((opredictedRUL, list(reversed(np.arange(xpredictedRUL[i])))))
opredictedRUL
The following figure shows the ZigZag pattern (first 1000 rows) for predicted and true RUL, and should show quite similar peak patterns:
mx = 1000
fig = plt.figure(figsize = (12, 8))
fig.add_subplot(1,2,1)
plt.plot(opredictedRUL[0:mx], color='blue')
plt.legend(['Predicted RUL'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylim(0, opredictedRUL[0:mx].max()+10)
plt.ylabel('RUL (cycles)')
fig.add_subplot(1,2,2)
plt.plot(otrueRUL[0:mx], color='purple')
plt.legend(['True RUL'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.ylabel('RUL (cycles)')
plt.ylim(0,otrueRUL[0:mx].max()+10)
plt.show()
16. Comparison predicted and true RUL
We can compare the predicted with the true RUL values using the following lineplot. Actually a lineplot is strictly spoken not valid here, since there are no measures inbetween the units. But it is for visualization purposes: the eye catches quickly that the predicted RULs are often a bit higher than the true values.
plt.figure(figsize = (16, 8))
plt.plot(RUL)
plt.plot(predictedRUL)
plt.xlabel('# Unit', fontsize=16)
plt.xticks(fontsize=16)
plt.ylabel('RUL', fontsize=16)
plt.yticks(fontsize=16)
plt.legend(['True RUL','Predicted RUL'], bbox_to_anchor=(0., 1.02, 1., .102), loc=3, mode="expand", borderaxespad=0)
plt.show()
Let us inspect the differences in a DataFrame. First we concatenate the true and predicted RUL:
df1 = pd.concat([pd.Series(RUL[0]), pd.Series(xpredictedRUL)], axis=1)
df1.columns = ['true', 'predicted']
And compute the difference score, which will show us if more values are positive or negative:
df1['diff'] = df1['predicted']-df1['true']
xpredictedRUL[0:5]
df1.head()
There is a tendency for the RUL testvalues being overestimated, as we can see in the histogram:
plt.hist(df1['diff'], bins=26, color="pink", edgecolor='brown', linewidth=1.2)
plt.axvline(0, color="red", linestyle='dashed', linewidth=1.6)
plt.show()
And the below table confirms this as well:
pd.DataFrame({'Count': [(df1['diff']<0).sum(), (df1['diff']==0).sum(), (df1['diff']>0).sum()]}, columns=['Count'], index=['Smaller', 'Zero', 'Larger'])
In many data science projects, slight overestimation is not a real problem. In the current case this can however be a risky business, with regards to the balance between maintenance costs and safety of aircrafts.
17. MSE
Finally, we may like to express the obtained result in a performance metric. MSE is often used in regression problems.
from sklearn.metrics import mean_squared_error
mse = mean_squared_error(RUL, xpredictedRUL)
print mse
Closing words
So far the results are quite reasonable for a first model in Keras, not denying that a lot can still be improved. Comparing with previous data, the present results are in a MSE range reported by others6. Several factors can still improve the result:
- Parameter tuning was not used and may give better results
- Cross-validation may control for overfitting
- Features may be improved by for example filtering out the noise
- It should be said that the chosen model may not be the best, other models may take into account the time-related aspects in the data (for example LSTM).
In addition, in a given project the performance metrics (MSE) may not be the main goal, insight into features that explain engine failure may be as important. That goal may be reached using algorithms that can give insight into feature importance such as random forest.
See you next time!
Hope you liked this blog. If you have any questions, suggestions, or see any errors, please contact me: rrighart@googlemail.com.
Notes & References
A. Saxena and K. Goebel (2008). "Turbofan Engine Degradation Simulation Data Set", NASA Ames Prognostics Data Repository (http://ti.arc.nasa.gov/project/prognostic-data-repository), NASA Ames Research Center, Moffett Field, CA
Windturbines. https://www.greentechmedia.com/articles/read/big-data-is-boosting-power-production-reducing-downtime-across-wind-fleets
A. Qin et al. (2017). Remaining useful life prediction for rotating machinery based on optimal degradation indicator. https://www.hindawi.com/journals/sv/2017/6754968/
Download a zip file. https://stackoverflow.com/questions/9419162/python-download-returned-zip-file-from-url
The name ZigZag should not be confused with the package "ZigZag" that is used to compute peaks in time-series. https://github.com/jbn/ZigZag/blob/master/zigzag_demo.ipynb
Jain, Kundu, Kumar Lad (2014). Prediction of remaining useful life of an aircraft engine under unknown initial wear. AIMTDR 2014.
(c) 2018, R. Righart | Website: https://rrighart.github.io | Email: rrighart@googlemail.com
Acknowledgements for CSS style go to jckantor