Information about the problem is at http://ti.arc.nasa.gov/tech/dash/pcoe/prognostic-data-repository/publications/#turbofan and original data is at http://ti.arc.nasa.gov/tech/dash/pcoe/prognostic-data-repository/#turbofan
The data was originally generated using the Commercial Modular Aero-Propulsion System Simulations (C-MAPPS) system.
The approach used in the turbofan engine degradation dataset was then used in the PHM08 challenge. Information about other research on the C-MAPSS data is available at https://www.phmsociety.org/sites/phmsociety.org/files/phm_submission/2014/phmc_14_063.pdf
In [1]:
import sys
import h2o
from h2o.estimators.glm import H2OGeneralizedLinearEstimator
from h2o.estimators.gbm import H2OGradientBoostingEstimator
from h2o.utils.shared_utils import _locate
import numpy as np
import pandas as pd
import seaborn as sns
import pykalman as pyk
sns.set()
doGridSearch = True
doKalmanSmoothing = False #unrelated to h2o, set true for demo
In [2]:
# Input files don't have column names
dependent_vars = ['RemainingUsefulLife']
index_columns_names = ["UnitNumber","Cycle"]
operational_settings_columns_names = ["OpSet"+str(i) for i in range(1,4)]
sensor_measure_columns_names =["SensorMeasure"+str(i) for i in range(1,22)]
input_file_column_names = index_columns_names + operational_settings_columns_names + sensor_measure_columns_names
# And we are going to add these columns
kalman_smoothed_mean_columns_names =["SensorMeasureKalmanMean"+str(i) for i in range(1,22)]
In [3]:
train = pd.read_csv("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/train_FD001.txt", sep=r"\s*", header=None,
names=input_file_column_names, engine='python')
test = pd.read_csv("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/test_FD001.txt", sep=r"\s*", header=None,
names=input_file_column_names, engine='python')
test_rul = pd.read_csv("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/RUL_FD001.txt", header=None, names=['RemainingUsefulLife'])
test_rul.index += 1 # set the index to be the unit number in the test data set
test_rul.index.name = "UnitNumber"
In [4]:
# Calculate the remaining useful life for each training sample based on last measurement being zero remaining
grouped_train = train.groupby('UnitNumber', as_index=False)
useful_life_train = grouped_train.agg({'Cycle' : np.max })
useful_life_train.rename(columns={'Cycle': 'UsefulLife'}, inplace=True)
train_wfeatures = pd.merge(train, useful_life_train, on="UnitNumber")
train_wfeatures["RemainingUsefulLife"] = -(train_wfeatures.UsefulLife - train_wfeatures.Cycle)
train_wfeatures.drop('UsefulLife', axis=1, inplace=True)
grouped_test = test.groupby('UnitNumber', as_index=False)
useful_life_test = grouped_test.agg({'Cycle' : np.max })
useful_life_test.rename(columns={'Cycle': 'UsefulLife'}, inplace=True)
test_wfeatures = pd.merge(test, useful_life_test, on="UnitNumber")
test_wfeatures["RemainingUsefulLife"] = -(test_wfeatures.UsefulLife - test_wfeatures.Cycle)
test_wfeatures.drop('UsefulLife', axis=1, inplace=True)
In [5]:
sns.set_context("notebook", font_scale=1.5)
p = sns.pairplot(train_wfeatures.query('UnitNumber < 10'),
vars=["RemainingUsefulLife", "SensorMeasure4", "SensorMeasure3",
"SensorMeasure9", "SensorMeasure8", "SensorMeasure13"], size=10,
hue="UnitNumber", palette=sns.color_palette("husl", 9));
sns.plt.show()
In [6]:
kalman_smoothed_mean_columns_names =["SensorMeasureKalmanMean"+str(i) for i in range(1,22)]
def calcSmooth(measures):
kf = pyk.KalmanFilter(initial_state_mean=measures[0], n_dim_obs=measures.shape[1])
(smoothed_state_means, smoothed_state_covariances) = kf.em(measures).smooth(measures)
return smoothed_state_means
def filterEachUnit(df):
dfout = df.copy()
for newcol in kalman_smoothed_mean_columns_names:
dfout[newcol] = np.nan
for unit in dfout.UnitNumber.unique():
sys.stdout.write('\rProcessing Unit: %d' % unit)
sys.stdout.flush()
unitmeasures = dfout[dfout.UnitNumber == unit][sensor_measure_columns_names]
smoothed_state_means = calcSmooth( np.asarray( unitmeasures ) )
dfout.loc[dfout.UnitNumber == unit, kalman_smoothed_mean_columns_names] = smoothed_state_means
sys.stdout.write('\rProcessing Unit: %d' % unit)
sys.stdout.flush()
sys.stdout.write('\rFinished\n')
sys.stdout.flush()
return dfout
In [7]:
# Get picky about the order of output columns
test_output_cols = index_columns_names + operational_settings_columns_names + sensor_measure_columns_names + \
kalman_smoothed_mean_columns_names
train_output_cols = test_output_cols + dependent_vars
if doKalmanSmoothing:
train_wkalman = filterEachUnit(train_wfeatures)
test_wkalman = filterEachUnit(test_wfeatures)
train_output = train_wkalman[train_output_cols]
test_output = test_wkalman[test_output_cols]
In [8]:
# Output the files, so we don't have to do the preprocessing again.
if doKalmanSmoothing:
train_output.to_csv("train_FD001_preprocessed.csv", index=False)
test_output.to_csv("test_FD001_preprocessed.csv", index=False)
test_rul.to_csv("rul_FD001_preprocessed.csv", index=True)
In [9]:
h2o.init()
In [10]:
#Pull Kalman-smoothed data if generated locally, or source from AWS
if doKalmanSmoothing:
train_hex = h2o.import_file(_locate("train_FD001_preprocessed.csv"))
test_hex = h2o.import_file(_locate("test_FD001_preprocessed.csv"))
else:
train_hex = h2o.import_file("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/train_FD001_preprocessed.csv")
test_hex = h2o.import_file("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/test_FD001_preprocessed.csv")
Use the operational settings and Kalman smoothed mean states as the independent features
Setup a fold column to great cross validation models from 90 units and cross validating on 10 units. This creates a 10-fold cross validation. The cross validation models are then used to create an ensemble model for predictions
In [11]:
xCols= operational_settings_columns_names + kalman_smoothed_mean_columns_names
yCol = dependent_vars
foldCol = "UnitNumberMod10"
train_hex[foldCol] = train_hex["UnitNumber"] % 10
In [12]:
def trainGLM(x, y, fold_column, training_frame, alpha=0.5, penalty=1e-5):
model = H2OGeneralizedLinearEstimator(family = "gaussian",alpha = [alpha], Lambda = [penalty])
model.train(x=x, y=y, training_frame=training_frame, fold_column=fold_column)
return model
def gridSearchGLM(x, y, fold_column, training_frame, alphas = [0,0.5,1], penalties=np.logspace(-3,0,num=4)):
results = []
for alpha in alphas:
for penalty in penalties:
results.append( trainGLM(x, y, fold_column, training_frame, alpha, penalty) )
return results
if doGridSearch:
glmModels = gridSearchGLM(xCols, yCol, foldCol, train_hex)
else:
# this is used to speed up the demonstration by just building the single model previously found
glmModels = [ trainGLM(xCols, yCol, foldCol, train_hex, alpha=1, penalty=0.01 )]
In [13]:
def extractBestModel(models):
bestMse = models[0].mse(xval=True)
result = models[0]
for model in models:
if model.mse(xval=True) < bestMse:
bestMse = model.mse(xval=True)
result = model
return result
bestModel = extractBestModel(glmModels)
bestModel
Out[13]:
In [14]:
def trainGBM(x, y, fold_column, training_frame, learning_rate=0.1, ntrees=50, max_depth=5):
model = H2OGradientBoostingEstimator(distribution = "gaussian",
learn_rate=learning_rate, ntrees=ntrees, max_depth=max_depth)
model.train(x=x, y=y, training_frame=training_frame, fold_column=fold_column)
return model
def gridSearchGBM(x, y, fold_column, training_frame, learning_rates = [0.1,0.03,0.01], ntrees=[10,30,100,300], max_depth=[1,3,5]):
results = []
for learning_rate in learning_rates:
for ntree in ntrees:
for depth in max_depth:
print "GBM: {learning rate: "+str(learning_rate)+"},{ntrees: "+str(ntree)+"},{max_depth: "+str(depth)+"}"
results.append( trainGBM(x, y, fold_column, training_frame, learning_rate=learning_rate, ntrees=ntree, max_depth=depth) )
return results
In [15]:
if doGridSearch:
#bmModels = gridSearchGBM(xCols, yCol, foldCol, train_hex,\
# learning_rates=[0.03,0.01,0.003], ntrees=[100,300,500], max_depth=[1,3,5])
#run the below line for fast demo
gbmModels = gridSearchGBM(xCols, yCol, foldCol, train_hex, learning_rates=[0.03,0.01], ntrees=[50,200], max_depth=[2,5])
else:
gbmModels = [trainGBM(xCols, yCol, foldCol, train_hex, \
ntrees=300, max_depth=5)]
In [16]:
bestGbmModel = extractBestModel(gbmModels)
Best model had depth 5, learning rate 0.01, and 300 trees
In [17]:
bestGbmModel.params
Out[17]:
Best GBM Model reported MSE on cross validation data as 1687, an improvement from GLM of 1954.
In [18]:
bestGbmModel
Out[18]:
In [19]:
train_hex["weights"] = 1
allModels = bestGbmModel.xvals
pred = sum([model.predict(train_hex) for model in allModels]) / len(allModels)
In [20]:
pred["actual"] = train_hex["RemainingUsefulLife"]
pred["unit"] = train_hex["UnitNumber"]
Ideally all points would be on the diagonal, indication prediction from data matched exactly the actual.
Also, it is important that the prediction gets more accurate the closer it gets to no useful life remaining.
Looking at a sample of the first 12 units.
Moved predictions from H2O to Python Pandas for plotting using Seaborn.
In [21]:
scored_df = pred.as_data_frame(use_pandas=True)
In [22]:
sns.set_context("notebook", font_scale=3)
g=sns.lmplot(x="actual",y="predict",hue="unit",col="unit",data=scored_df[scored_df.unit < 13],col_wrap=3,fit_reg=False, size=10)
ticks = np.linspace(-300,100, 5)
g = (g.set_axis_labels("Remaining Useful Life", "Predicted Useful Life")
.set(xlim=(-325, 125), ylim=(-325, 125),
xticks=ticks, yticks=ticks))
In [23]:
np.linspace(-300,100, 5)
Out[23]:
In [24]:
testPreds = sum([model.predict(test_hex) for model in allModels]) / len(allModels)
Append the original index information (Cycle and UnitNumber) to the predicted values so we have them later.
In [25]:
testPreds["Cycle"] = test_hex["Cycle"]
testPreds["UnitNumber"] = test_hex["UnitNumber"]
Move the predictions over to Python Pandas for final analysis and scoring
In [26]:
testPreds_df = testPreds.as_data_frame(use_pandas=True)
Load up the actual Remaining Useful Life information.
In [27]:
if doKalmanSmoothing:
actual_RUL = pd.read_csv(_locate("rul_FD001_preprocessed.csv"))
else:
actual_RUL = pd.read_csv("http://h2o-public-test-data.s3.amazonaws.com/bigdata/laptop/CMAPSSData/rul_FD001_preprocessed.csv")
The final scoring used in the competition is based on a single value per unit. We extract the last three predictions and use the mean of those (simple aggregation) and put the prediction back from remaining useful life in T-minus format to cycles remaining (positive).
In [28]:
def aggfunc(x):
return np.mean( x.order().tail(3) )
grouped_by_unit_preds = testPreds_df.groupby("UnitNumber", as_index=False)
predictedRUL = grouped_by_unit_preds.agg({'predict' : aggfunc })
predictedRUL.predict = -predictedRUL.predict
Add the prediction to the actual data frame, and use the scoring used in the PHMO8 competition (more penality for predicting more useful life than there is actual).
In [29]:
final = pd.concat([actual_RUL, predictedRUL.predict], axis=1)
In [30]:
def rowScore(row):
d = row.predict-row.RemainingUsefulLife
return np.exp( -d/10 )-1 if d < 0 else np.exp(d/13)-1
rowScores = final.apply(rowScore, axis=1)
This is the final score using PHM08 method of scoring.
In [31]:
sum(rowScores)
Out[31]:
In [32]:
sns.set_context("notebook", font_scale=1.25)
sns.regplot("RemainingUsefulLife", "predict", data=final, fit_reg=False);