Прогнозирование простоев (профилактическое обслуживание) оборудования с использованием регрессионных моделей для прогнозирования простоя…

Название может показаться слишком длинным, но проект - это реальный вариант использования, с которым я столкнулся во время своего опыта. Большинство производственных компаний изменили свою стратегию получения прибыли. Они продают свою продукцию по себестоимости в обмен на то, что хотели получать доход и прибыль, продавая свои услуги и техническое обслуживание, с которыми машины сталкиваются в течение своей жизни.

Чтобы выполнить финансовый план после выхода на рынок после выхода на рынок, всем этим компаниям необходимо предсказать, когда именно оборудование выйдет из строя. До появления достижений в области машинного обучения большинство производственных компаний используют для обзвона своих клиентов через определенные месяцы. с даты покупки машины или по фиксированной шкале времени с даты последнего технического обслуживания. В этом случае клиент может подумать о вариантах или местных точках обслуживания (которые считаются поддельными).

Чтобы удовлетворить вышеперечисленные требования, аналитические группы, работающие над созданием систем раннего оповещения, для прогнозирования условий простоя до того, как он действительно произойдет в оборудовании. Применяя подход, основанный на анализе данных, компании-производители могут обращаться к своим клиентам за месяц до того, как оборудование достигает срока обслуживания, и таким образом превращать их в своих послепродажных клиентов, не теряя доли рынка в бизнесе.

Это была одна из самых сложных проблем, с которыми я имел дело. Нам были предоставлены данные 45 датчиков, которые собирают различные показатели оборудования, такие как температура, значения давления, состояние аварийного напряжения, реактивная мощность в агрегате, угловое расположение валов, ротор и статор, максимальный крутящий момент и т. Д. Зависимой переменной является количество дней Оставшееся до того, как оборудование вышло из строя из-за «регресса», а количество дней, оставшихся на техническое обслуживание, меньше и больше 30 дней для классификации.

Теперь теория готова, давайте перейдем к решению проблем с использованием традиционных алгоритмов машинного обучения и LSTM для классификации.

Загрузка всех необходимых пакетов

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import keras
from sklearn.preprocessing import MinMaxScaler
from sklearn.metrics import confusion_matrix,accuracy_score
from keras.models import Sequential
from keras.layers import Dense, Dropout, LSTM, Activation
from keras.callbacks import EarlyStopping
import matplotlib.pyplot as plt
plt.style.use(‘ggplot’)
%matplotlib inline

Загрузка данных в предпочтительную среду в примере, который я показываю, загружается в python3.6.

train=pd.read_csv(‘../train.txt’,sep=’ ‘,header=None),axis=1)
col_names = [‘UnitNumber’, ‘Cycle’, ‘Op_Setting_1’, ‘Op_Setting_2’, ‘Op_Setting_3’, ‘Sensor_1’, ‘Sensor_2’, ‘Sensor_3’, ‘Sensor_4’, ‘Sensor_5’, ‘Sensor_6’, ‘Sensor_7’, ‘Sensor_8’, ‘Sensor_9’, ‘Sensor_10’, ‘Sensor_11’, ‘Sensor_12’, ‘Sensor_13’, ‘Sensor_14’, ‘Sensor_15’, ‘Sensor_16’, ‘Sensor_17’, ‘Sensor_18’, ‘Sensor_19’, ‘Sensor_20’, ‘Sensor_21’]
train.columns=col_names
print(‘Shape of Train dataset: ‘,dataset_train.shape)
train.head()

Код возвращает заголовок фрейма данных с 5 строками и назначенными столбцами.

Shape of Train data set:  (20631, 26)

test=pd.read_csv(‘../test.txt’,sep=’‘,header=None).drop([26,27],axis=1)
test.columns=col_names
print(‘Shape of Test dataset: ‘,dataset_train.shape)
test.head()

Код возвращает заголовок фрейма данных с 5 строками и назначенными столбцами.

Shape of Test dataset:  (20631, 26)

df=pd.concat([train,test])
max_cycle = df.groupby(‘UnitNumber’)[‘Cycle’].max().reset_index()
max_cycle.columns = [‘UnitNumber’, ‘MaxOfCycle’]

df_merged = df.merge(max_cycle, left_on=’UnitNumber’, right_on=’UnitNumber’, how=’inner’)

Target_Remaining_Useful_Life = df_merged[“MaxOfCycle”] — df_merged[“Cycle”]
df_with_target = df_merged[“Target_Remaining_Useful_Life”] = Target_Remaining_Useful_Life

df_with_target = df_merged.drop(“MaxOfCycle”, axis=1)
df_with_target[df_with_target[‘UnitNumber’] == 1].head(5)
target_var = [‘Target_Remaining_Useful_Life’]
index_columns_names = [“UnitNumber”,”Cycle”]
op_settings_columns = [“Op_Setting_”+str(i) for i in range(1,4)]
sensor_columns =[“Sensor_”+str(i) for i in range(1,22)]
print(df_with_target.shape)
leakage_to_drop = ['UnitNumber', 'Cycle', 'Op_Setting_1', 'Op_Setting_2', 'Op_Setting_3','Sensor_18','Sensor_19']  
df_no_leakage = df_with_target.drop(leakage_to_drop, axis = 1)
print(df_no_leakage.shape)
# set up features and target variable 
y = df_no_leakage['Target_Remaining_Useful_Life']
X = df_no_leakage.drop(['Target_Remaining_Useful_Life'], axis = 1)
df=df_no_leakage
corr=df.corr()
plt.rcParams[“font.size”] = 12
plt.rcParams[“figure.figsize”] = (len(corr),len(corr))
plt.subplot(1, 1, 1)
sns.heatmap(corr, square=True, vmax=1, vmin=-1, center=0,annot=True,fmt=”1.1f”)
plt.show

Использование техники ранжирования характеристик для определения лучших характеристик

from sklearn import ensemble
rf = ensemble.RandomForestRegressor()
single_rf = ensemble.RandomForestRegressor(n_estimators = 200, max_depth = 15)
single_rf.fit(X, y)
y_pred = single_rf.predict(X)
print(“complete”)
import matplotlib.pyplot as plt
importances = single_rf.feature_importances_
indices = np.argsort(importances)[::-1]
feature_names = X.columns 
f, ax = plt.subplots(figsize=(11, 9))
plt.title(“Feature ranking”, fontsize = 20)
plt.bar(range(X.shape[1]), importances[indices], color=”b”, align=”center”)
plt.xticks(range(X.shape[1]), indices) #feature_names, rotation=’vertical’)
plt.xlim([-1, X.shape[1]])
plt.ylabel(“importance”, fontsize = 18)
plt.xlabel(“index of the feature”, fontsize = 18)
plt.show()
# list feature importance
important_features = pd.Series(data=single_rf.feature_importances_,index=X.columns)
important_features.sort_values(ascending=False,inplace=True)
print(important_features.head(10))

print(train_no_leakage.shape)
vars_to_drop = [“Sensor_”+str(i) for i in [5, 15, 9, 17, 4]]
df_final = df_no_leakage.drop(vars_to_drop, axis = 1)
print(df_final.shape)
from sklearn import preprocessing
categorical = df_final.select_dtypes(include=['object'])
numeric = df_final.select_dtypes(exclude=['object'])
print(categorical.columns.values)

for name, values in categorical.items():
    print(name)
    dummies = pd.get_dummies(values.str.strip(), prefix = name, dummy_na=True)
    numeric = pd.concat([numeric, dummies], axis=1)

for name in numeric:
    print(name)
    if pd.isnull(numeric[name]).sum() > 0:
        numeric["%s_mi" % (name)] = pd.isnull(numeric[name])
        median = numeric[name].median()
        numeric[name] = numeric[name].apply(lambda x: median if pd.isnull(x) else x)
y = numeric['Target_Remaining_Useful_Life']
X = numeric.drop(['Target_Remaining_Useful_Life'], axis = 1)

Применение регрессии случайного леса для прогнозирования оставшегося срока полезного использования до технического обслуживания

import numpy as np
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1234)
# choose the model
from sklearn.ensemble import RandomForestRegressor
rf = ensemble.RandomForestRegressor()
# set up 5-fold cross-validation
from sklearn import model_selection
cv = model_selection.KFold(5)
# pipeline standardization and model
from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps=[(‘standardize’, preprocessing.StandardScaler())
 , (‘model’, rf) ])
# tune the model
my_min_samples_leaf = [2, 10, 25, 50, 100]
my_max_depth = [7, 8, 9, 10, 11, 12]
# run the model using gridsearch, select the model with best search
from sklearn.model_selection import GridSearchCV
optimized_rf = GridSearchCV(estimator=pipeline
 , cv=cv
 , param_grid =dict(model__min_samples_leaf = my_min_samples_leaf, model__max_depth = my_max_depth)
 , scoring = ‘neg_mean_squared_error’
 , verbose = 1
 , n_jobs = -1
 )
optimized_rf.fit(X_train, y_train)
# show the best model estimators
print(optimized_rf.best_estimator_)
# evaluate metrics on holdout
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
y_pred = optimized_rf.predict(X_test)
print(“Random Forest Mean Squared Error: “, mean_squared_error(y_test, y_pred))
print(“Random Forest Mean Absolute Error: “, mean_absolute_error(y_test, y_pred))
print(“Random Forest r-squared: “, r2_score(y_test, y_pred))

Применение ElasticNet, которое учитывает регрессию Лассо и Ридж при уменьшении параметров

import numpy as np
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1234)
# choose the model
from sklearn.linear_model import ElasticNet
glm_net = ElasticNet()
# set up 5-fold cross-validation
from sklearn import model_selection
cv = model_selection.KFold(5)
# pipeline standardization and model
from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps=[(‘standardize’, preprocessing.StandardScaler())
 , (‘model’, glm_net) ])
# tune the model
my_alpha = np.linspace(.01, 1, num=5)
my_l1_ratio = np.linspace(.01, 1, num=3)
# run the model using gridsearch, select the model with best search
from sklearn.model_selection import GridSearchCV
optimized_glm_net = GridSearchCV(estimator=pipeline
 , cv=cv
 , param_grid =dict(model__l1_ratio = my_l1_ratio, model__alpha = my_alpha)
 , scoring = ‘neg_mean_squared_error’
 , verbose = 1
 , n_jobs = -1
 )
optimized_glm_net.fit(X_train, y_train)
# show the best model estimators
print(optimized_glm_net.best_estimator_)
# evaluate metrics on holdout
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
y_pred = optimized_glm_net.predict(X_test)
print(“GLM Elastic Net Mean Squared Error: “, mean_squared_error(y_test, y_pred))
print(“GLM Elastic Net Mean Absolute Error: “, mean_absolute_error(y_test, y_pred))
print(“GLM Elastic Net r-squared: “, r2_score(y_test, y_pred))

Применение регрессора опорных векторов

import numpy as np
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1234)
# choose the model
from sklearn import svm
from sklearn.svm import SVR
svm = svm.SVR()
# set up 5-fold cross-validation
from sklearn import model_selection
cv = model_selection.KFold(5)
# pipeline standardization and model
from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps=[(‘standardize’, preprocessing.StandardScaler())
 , (‘model’, svm) ])
# tune the model
my_C = [1]
my_epsilon = [.05, .1, .15]
# run the model using gridsearch, select the model with best search
from sklearn.model_selection import GridSearchCV
optimized_svm = GridSearchCV(estimator=pipeline
 , cv=cv
 , param_grid =dict(model__C = my_C, model__epsilon = my_epsilon)
 , scoring = ‘neg_mean_squared_error’
 , verbose = 1
 , n_jobs = -1
 )
optimized_svm.fit(X_train, y_train)
# show the best model estimators
print(optimized_svm.best_estimator_)
# evaluate metrics on holdout
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
y_pred = optimized_svm.predict(X_test)
print(“SVM Mean Squared Error: “, mean_squared_error(y_test, y_pred))
print(“SVM Mean Absolute Error: “, mean_absolute_error(y_test, y_pred))
print(“SVM r-squared: “, r2_score(y_test, y_pred))

Применение алгоритма повышения градиента

import numpy as np
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1234)
# choose the model
from sklearn.ensemble import GradientBoostingRegressor
gb = ensemble.GradientBoostingRegressor()
# set up 5-fold cross-validation
from sklearn import model_selection
cv = model_selection.KFold(5)
# pipeline standardization and model
from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps=[(‘standardize’, preprocessing.StandardScaler())
 , (‘model’, gb) ])
# tune the model
my_alpha = [.5, .75, .9]
my_n_estimators= [500]
my_learning_rate = [0.005, .01]
my_max_depth = [4, 5, 6]
# run the model using gridsearch, select the model with best search
from sklearn.model_selection import GridSearchCV
optimized_gb = GridSearchCV(estimator=pipeline
 , cv=cv
 , param_grid =dict(model__max_depth = my_max_depth, model__n_estimators = my_n_estimators,
 model__learning_rate = my_learning_rate, model__alpha = my_alpha)
 , scoring = ‘neg_mean_squared_error’
 , verbose = 1
 , n_jobs = -1
 )
optimized_gb.fit(X_train, y_train)
# show the best model estimators
print(optimized_gb.best_estimator_)
# evaluate metrics on holdout
from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score
y_pred = optimized_gb.predict(X_test)
print(“Gradient Boosting Mean Squared Error: “, mean_squared_error(y_test, y_pred))
print(“Gradient Boosting Mean Absolute Error: “, mean_absolute_error(y_test, y_pred))
print(“Gradient Boosting r-squared: “, r2_score(y_test, y_pred))

fig, ax = plt.subplots()
ax.scatter(y_test, y_pred, edgecolors=(0, 0, 0))
ax.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], ‘k — ‘, lw=4)
ax.set_xlabel(‘Actual RUL’)
ax.set_ylabel(‘Predicted RUL’)
ax.set_title(‘Remaining Useful Life Actual vs. Predicted’)
plt.show()

cycles = 15
df_no_leakage[‘Target_15_Cycles’] = np.where(df_no_leakage[‘Target_Remaining_Useful_Life’] <= cycles, 1, 0 )
df_no_leakage.tail(5)

print(train_no_leakage.shape)
vars_to_drop = [“Sensor_”+str(i) for i in [5, 15, 9, 17, 4]]
target_to_drop = [‘Target_Remaining_Useful_Life’]
train_final = train_no_leakage.drop(vars_to_drop, axis = 1)
train_final = train_no_leakage.drop(target_to_drop, axis = 1)
train_final.tail()

from sklearn import preprocessing
categorical = train_final.select_dtypes(include=[‘object’])
numeric = train_final.select_dtypes(exclude=[‘object’])
print(categorical.columns.values)
# create dummy variables (if any categorical fields)
for name, values in categorical.items():
 print(name)
 dummies = pd.get_dummies(values.str.strip(), prefix = name, dummy_na=True)
 numeric = pd.concat([numeric, dummies], axis=1)
# imputation (if any NULL values)
for name in numeric:
 print(name)
 if pd.isnull(numeric[name]).sum() > 0:
 numeric[“%s_mi” % (name)] = pd.isnull(numeric[name])
 median = numeric[name].median()
 numeric[name] = numeric[name].apply(lambda x: median if pd.isnull(x) else x)
y = numeric[‘Target_15_Cycles’]
X = numeric.drop([‘Target_15_Cycles’], axis = 1)

import numpy as np
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=1234)
# choose the model
from sklearn import ensemble
from sklearn.ensemble import RandomForestClassifier
rf = ensemble.RandomForestClassifier()
# set up 5-fold cross-validation
from sklearn import model_selection
cv = model_selection.KFold(5)
# pipeline standardization and model
from sklearn.pipeline import Pipeline
pipeline = Pipeline(steps=[(‘standardize’, preprocessing.StandardScaler())
 , (‘model’, rf) ])
# tune the model
my_min_samples_leaf = [2, 25, 50]
my_max_depth = [8, 9, 10, 12]
# run the model using gridsearch, select the model with best search
from sklearn.model_selection import GridSearchCV
optimized_rf = GridSearchCV(estimator=pipeline
 , cv=cv
 , param_grid =dict(model__min_samples_leaf = my_min_samples_leaf, model__max_depth = my_max_depth)
 , scoring = ‘roc_auc’
 , verbose = 1
 , n_jobs = -1
 )
optimized_rf.fit(X_train, y_train)
# show the best model estimators
y_pred_proba = optimized_rf.predict_proba(X_test)[:, 1]
y_pred = optimized_rf.predict(X_test)
print(optimized_rf.best_estimator_)

from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, confusion_matrix
print(“Confusion Matrix:”)
print(confusion_matrix(y_test, y_pred))
from sklearn.metrics import classification_report
print(“Random Forest Accuracy: “+”{:.1%}”.format(accuracy_score(y_test, y_pred)));
print(“Random Forest Precision: “+”{:.1%}”.format(precision_score(y_test, y_pred)));
print(“Random Forest Recall: “+”{:.1%}”.format(recall_score(y_test, y_pred)));
print(“Classification Report:”)
print(classification_report(y_test, y_pred))
from sklearn import metrics
fpr, tpr, threshold = metrics.roc_curve(y_test, y_pred_proba)
roc_auc = metrics.auc(fpr, tpr)
import matplotlib.pyplot as plt
plt.title(‘Receiver Operating Characteristic’)
plt.plot(fpr, tpr, ‘b’, label = ‘AUC = %0.2f’ % roc_auc)
plt.legend(loc = ‘lower right’)
plt.plot([0, 1], [0, 1],’r — ‘)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel(‘True Positive Rate’)
plt.xlabel(‘False Positive Rate’)
plt.show()

def gen_sequence(id_df, seq_length, seq_cols):
 df_zeros=pd.DataFrame(np.zeros((seq_length-1,id_df.shape[1])),columns=id_df.columns)
 id_df=df_zeros.append(id_df,ignore_index=True)
 data_array = id_df[seq_cols].values
 num_elements = data_array.shape[0]
 lstm_array=[]
 for start, stop in zip(range(0, num_elements-seq_length), range(seq_length, num_elements)):
 lstm_array.append(data_array[start:stop, :])
 return np.array(lstm_array)
# function to generate labels
def gen_label(id_df, seq_length, seq_cols,label):
 df_zeros=pd.DataFrame(np.zeros((seq_length-1,id_df.shape[1])),columns=id_df.columns)
 id_df=df_zeros.append(id_df,ignore_index=True)
 data_array = id_df[seq_cols].values
 num_elements = data_array.shape[0]
 y_label=[]
 for start, stop in zip(range(0, num_elements-seq_length), range(seq_length, num_elements)):
 y_label.append(id_df[label][stop])
 return np.array(y_label)
seq_length=50
seq_cols=numeric.columns.values
cycles = 15
df_with_target['Target_15_Cycles'] = np.where(df_with_target['Target_Remaining_Useful_Life'] <= cycles, 1, 0 )
df_with_target.tail(5)

sc=MinMaxScaler()
df_with_target[seq_cols]=sc.fit_transform(df_with_target[seq_cols])
X=np.concatenate(list(list(gen_sequence(df_with_target[df_with_target['UnitNumber']==id], seq_length, seq_cols)) for id in df_with_target['UnitNumber'].unique()))
print(X.shape)
# generate y_train
y=np.concatenate(list(list(gen_label(df_with_target[df_with_target['UnitNumber']==id], 50, seq_cols,'Target_15_Cycles')) for id in df_with_target['UnitNumber'].unique()))
print(y.shape)
X_train, X_test, y_train, y_test = train_test_split(X, y , test_size=0.3, random_state=0)

Применение LSTM с помощью Kera’s

nb_features =X_train.shape[2]
timestamp=seq_length
model = Sequential()
model.add(LSTM(
 input_shape=(timestamp, nb_features),
 units=100,
 return_sequences=True))
model.add(Dropout(0.2))
model.add(LSTM(
 units=50,
 return_sequences=False))
model.add(Dropout(0.2))
model.add(Dense(units=1, activation=’sigmoid’))
model.compile(loss=’binary_crossentropy’, optimizer=’adam’, metrics=[‘accuracy’])
model.summary()

model.fit(X_train, y_train, epochs=10, batch_size=200, validation_split=0.05, verbose=1,
 callbacks = [EarlyStopping(monitor=’val_loss’, min_delta=0, patience=0, verbose=0, mode=’auto’)])

scores = model.evaluate(X_train, y_train, verbose=1, batch_size=200)
print(‘Accurracy: {}’.format(scores[1]))

y_pred=model.predict_classes(X_test)
y_pred_proba = model.predict_proba(X_test)
print(‘Accuracy of model on test data: ‘,accuracy_score(y_test,y_pred))
print(‘Confusion Matrix: \n’,confusion_matrix(y_test,y_pred))

from sklearn.metrics import roc_auc_score, accuracy_score, precision_score, recall_score, confusion_matrix
print(“Confusion Matrix:”)
print(confusion_matrix(y_test, y_pred))
from sklearn.metrics import classification_report
print(“LSTM Accuracy: “+”{:.1%}”.format(accuracy_score(y_test, y_pred)));
print(“LSTM Precision: “+”{:.1%}”.format(precision_score(y_test, y_pred)));
print(“LSTM Recall: “+”{:.1%}”.format(recall_score(y_test, y_pred)));
print(“Classification Report:”)
print(classification_report(y_test, y_pred))
from sklearn import metrics
fpr, tpr, threshold = metrics.roc_curve(y_test, y_pred_proba)
roc_auc = metrics.auc(fpr, tpr)
import matplotlib.pyplot as plt
plt.title(‘Receiver Operating Characteristic’)
plt.plot(fpr, tpr, ‘b’, label = ‘AUC = %0.2f’ % roc_auc)
plt.legend(loc = ‘lower right’)
plt.plot([0, 1], [0, 1],’r — ‘)
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel(‘True Positive Rate’)
plt.xlabel(‘False Positive Rate’)
plt.show()

С помощью подхода, которого мы придерживались, мы смогли продемонстрировать, что ранние сигналы могут быть обнаружены, поэтому компании могут поддерживать достаточный запас запасов и выделять нужных специалистов по обслуживанию для удовлетворения и удовлетворения потребностей в обслуживании и техническом обслуживании.

Я надеюсь, что эта статья поможет большинству команд по всему миру справиться с профилактическим обслуживанием.