Почти все области жизни генерируют данные с огромной скоростью. Затем эти данные используются для принятия прибыльных/безопасных и эффективных решений. Стандартный производственный процесс не является исключением. Данные, генерируемые различными станками во время обработки, могут использоваться для различных задач, таких как профилактическое обслуживание, отслеживание производительности станков, прогнозирование результатов и т. д. Одним из наиболее фундаментальных методов, используемых во многих из этих приложений, является прогнозирование временных рядов и обнаружение аномалий. Модели прогнозирования позволяют предсказать, как данные будут выглядеть в будущем. Алгоритм обнаружения аномалий сообщает нам, есть ли в данных какие-либо точки, которые являются ненормальными или аномальными в наборе данных. Эти данные могут быть очень важными при принятии различных решений. Чтобы продемонстрировать эти методы, я использую Набор данных износа при фрезеровании НАСА [1]. Весь код для анализа можно получить на GitHub.
О наборе данных
Эксперимент состоит в многократном фрезеровании блока фиксированных размеров до тех пор, пока используемый инструмент не изнашивается. Всего имеется 8 различных тестовых случаев, каждый из которых повторяется дважды. Различные тестовые примеры представляют собой комбинации различной глубины резания (1,5 мм или 0,75 мм), скоростей подачи (0,25 или 0,5) и, наконец, материала режущего инструмента (сталь или чугун).
Каждый случай генерирует следующие данные:
(a) VB-Величина износа задней поверхности.
(b) smcAC-Ток двигателя шпинделя переменного тока
(c) smcDC- ток двигателя шпинделя постоянного тока
(d) vib_table — данные о вибрации, собранные датчиком вибрации, удерживаемым в фиксированном положении стола.
(e) vib_spindle — данные о вибрации, собранные датчиком вибрации, хранящимся на шпинделе.
(f) AE_table — Акустическая эмиссия за столом
(g) AE_spindle — акустическая эмиссия на шпинделе
Дополнительную информацию можно получить из [1]
Исследовательский анализ данных
Давайте посмотрим, как различные параметры влияют на износ инструмента. График времени выполнения для различных экспериментов выглядит следующим образом:
Из графика видно, что вариант 11 имеет самое продолжительное время выполнения эксперимента. Это наблюдение может быть подкреплено существующими знаниями в области машиностроения.
Глубина резания в этом случае находится на нижнем пределе (0,75). Меньшая глубина резания означает меньшие силы резания, а значит, инструмент будет изнашиваться медленнее.
Значение подачи находится на нижнем конце (0,25). При более низких значениях подачи создаются меньшие динамические силы, что также снижает скорость износа инструмента.
Наконец, материал режущего инструмента корпуса 11 — чугун (1). Чугун способен поглощать больше ударов по сравнению со сталью из-за более высокого содержания углерода. Обладая большей способностью поглощать удары, режущий инструмент способен лучше поглощать динамические нагрузки и, как следствие, служит дольше.
С другой стороны, давайте рассмотрим случай 16, который представляет собой случай с наименьшим временем выполнения, который был успешно выполнен (случай с наименьшим временем выполнения, случай 6 игнорируется, поскольку он имеет только одну точку данных и не имеет статистического значения). .
Глубина резания выше (1,5), скорость подачи выше (0,5), а материал — сталь (2). В результате этот инструмент изнашивается намного быстрее, чем тот, который использовался в случае 11.
Поскольку случай 11 имеет максимальное количество точек данных, это наиболее подходящий случай для демонстрации прогнозирования и обнаружения аномалий.
Давайте посмотрим, как различные параметры ведут себя в зависимости от времени для этого конкретного эксперимента.
Поведение датчиков акустической эмиссии на шпинделе и столе в чем-то похоже, рассмотрим подробнее.
Из графиков можно сделать вывод, что динамика акустической эмиссии обоих датчиков более или менее одинакова. Датчик акустической эмиссии на шпинделе испытывает несколько большую акустическую эмиссию по величине. Это имеет смысл, так как он намного ближе к тому месту, где происходит фактическая обработка, по сравнению с датчиком стола, который расположен дальше.
Прогнозирование: инженерия данных и базовые модели
Попробуем спрогнозировать vib_spindle для эксперимента 11.
Можно заметить, что временной ряд неравномерен. Этот временной ряд необходимо сделать регулярным для дальнейшего анализа. Этот временной ряд делается регулярным путем его преобразования в формат даты и времени, а затем с помощью линейной интерполяции для заполнения промежуточных значений времени. Этот процесс преобразования нерегулярного временного ряда в обычный называется повторной выборкой. Это достигается следующими строками кода:
#Converting to date-time format
df3_vib_spindle.set_index(pd.to_datetime(df3_vib_spindle['time'], unit='s'), inplace=True)
df3_vib_spindle.drop('time', axis=1, inplace=True)
#Resampling
df3_vib_spindle_resampled = df3_vib_spindle.resample('1S').interpolate(method = 'linear')
Пересчитанный временной ряд теперь выглядит следующим образом:
Давайте проверим, есть ли в этом временном ряду сезонные или остаточные компоненты. Это можно наблюдать, выполняя аддитивную сезонную декомпозицию временного ряда.
Можно заметить, что нет сезонной или остаточной составляющей.
Прежде чем приступить к фактическому прогнозированию, важно иметь стандарт для оценки эффективности нашей модели прогнозирования. Этот стандарт устанавливается путем создания базовых моделей. Для этого набора данных были созданы две базовые модели:
(a) Среднее значение исторических данных
(b) Последнее значение исторических данных
Здесь я расскажу о случае (b). Вариант (а) можно найти в моем блокноте Jupyter Notebook.
Базовая модель: последнее значение исторических данных
Показателем оценки, выбранным для всех моделей, является «средняя ошибка в процентах» (MAPE). MAPE — популярный показатель, используемый в моделях прогнозирования, поскольку он выражает результаты в виде процентной ошибки и хорошо подходит для данных, которые почти не имеют нулевых значений.
В Python MAPE можно рассчитать с помощью следующей функции:
def mape(y_true, y_pred):
return np.mean((np.abs(y_true - y_pred)/y_true)) * 100
Для разработки и оценки модели модель делится на обучающие (исторические) данные и тестовый набор с разделением 80:20.
Код для создания базовой модели последнего значения приведен ниже:
#Storing the last value of the training set last = train['vib_spindle'].iloc[-1] #Adding 'prediction' to test set test.loc[:, 'pred_last'] = last
Затем рассчитывается MAPE для нашей базовой модели последнего значения:
mape_last = mape(test['vib_spindle'], test['pred_last']) mape_last
Расчет MAPE для базовой модели дает погрешность 1,77%.
Производительность базовой модели на данных тестового набора можно визуализировать следующим образом:
fig, ax = plt.subplots() ax.plot(train['vib_spindle'], 'g-', label = 'Train') ax.plot(test['vib_spindle'], 'b-', label = 'Test') ax.plot(test['pred_last'], 'r--', label = 'Predicted') ax.legend(loc = 2) plt.plot()
Теперь у нас есть базовая модель, с которой можно сравнить наши модели прогнозирования.
Прогнозирование: стационарность
Для большинства моделей статистического прогнозирования требуются стационарные временные ряды[2]. Стационарный временной ряд определяется как временной ряд, статистические свойства которого (среднее значение, дисперсия, автокорреляция и т. д.) не меняются с течением времени.
Обычной практикой является проверка стационарности временных рядов, которые необходимо прогнозировать. Если она не стационарна, то ее делают стационарной путем применения различных преобразований.
Давайте проверим наш временной ряд, чтобы увидеть, является ли он стационарным или нет. Одним из наиболее распространенных тестов на стационарность является расширенный тест Дикки Фуллера (ADF).
from statsmodels.tsa.stattools import adfuller
ADF_result = adfuller(df3_vib_spindle_resampled['vib_spindle'])
print(f'ADF statistic: {ADF_result[0]}' )
print(f'p-value: {ADF_result[1]}' )
ADF statistic: -1.9066453259812026
p-value: 0.32894142658345893
Поскольку статистика ADF не является большим отрицательным числом, а значение p больше 0,05, мы не можем отвергнуть нулевую гипотезу, а временной ряд не является стационарным.
Одно простое преобразование, позволяющее сделать временной ряд стационарным, — это разность первого порядка. При дифференцировании первого порядка текущее значение заменяется разницей между текущим значением и предыдущим значением.
#First order differencing df4 = df3_vib_spindle_resampled['vib_spindle'].diff() #Dropping first value after transformation df4 = df4[1:]
Давайте посмотрим, помогло ли преобразование сделать временные ряды стационарными.
ADF_result_diff = adfuller(df4['vib_spindle'])
print(f'ADF statistic: {ADF_result_diff[0]}' )
print(f'p-value: {ADF_result_diff[1]}' )
ADF statistic: -7.295304581462539
p-value: 1.3829646588434805e-10
Статистика ADF теперь представляет собой большое отрицательное число, а p-значение меньше 0,05. Следовательно, теперь мы можем отвергнуть нулевую гипотезу. Следовательно, временной ряд теперь стационарен и может использоваться для моделирования.
Прогнозирование: АКФ
Чтобы определить порядок теперь стационарных временных рядов, мы теперь построим график автокорреляционной функции (ACF). Автокорреляция измеряет линейную зависимость между запаздывающими значениями временного ряда.
from statsmodels.graphics.tsaplots import plot_acf plot_acf(df4['vib_spindle'], lags = 50)
После 1 значения запаздывающих признаков резко падают и находятся в пределах доверительного интервала. Следовательно, можно заметить, что временной ряд является временным рядом корреляции первого порядка. Основываясь на этих данных, сильным кандидатом для модели прогнозирования является модель скользящего среднего.
Прогнозирование: модель скользящего среднего
Модель скользящего среднего реализована с моделью SARIMAX через библиотеку statsmodel.
from statsmodels.tsa.statespace.sarimax import SARIMAX
def rolling_forecast(df:pd.DataFrame, train_len: int, horizon: int, window: int):
total_len = train_len + horizon
pred_MA = []
for i in range(train_len, total_len, window):
model = SARIMAX(df[:i], order = (0,0,1))
res = model.fit(disp = False)
predictions = res.get_prediction(0, i + window - 1)
oos_pred = predictions.predicted_mean.iloc[-window:]
pred_MA.extend(oos_pred)
return pred_MA
pred_df = test_diff.copy()
train_len = len(train_diff)
horizon = len(test_diff)
window = 1
pred_MA = rolling_forecast(df4['vib_spindle'], train_len, horizon, window)
pred_df['pred_MA'] = pred_MA
pred_df.head()
df3_vib_spindle_resampled['prediction_MA'] = pd.Series()
df3_vib_spindle_resampled['prediction_MA'] [84 :] = df3_vib_spindle_resampled['vib_spindle'].iloc[84] + pred_df['pred_MA'].cumsum()
MAPE для модели скользящего среднего получается следующим образом:
from sklearn.metrics import mean_absolute_percentage_error mape_MA_undiff = mean_absolute_percentage_error(df3_vib_spindle_resampled['vib_spindle'].iloc[84 :],df3_vib_spindle_resampled['prediction_MA'][84:]) print(mape_MA_undiff) 0.011796312171518794
MAPE для нашей модели скользящего среднего составляет 0,0118%, что намного лучше, чем наша базовая модель. Производительность модели скользящего среднего на тестовом наборе визуализирована ниже:
fig, ax = plt.subplots() ax.plot(df3_vib_spindle_resampled['vib_spindle'], 'b-', label = 'actual') ax.plot(df3_vib_spindle_resampled['prediction_MA'], 'k--', label = 'MA(1)') ax.legend()
Обнаружение аномалий
Обнаружение аномалий определяется как процесс изучения временного ряда и обнаружения точек данных, которые являются выбросами [4]. Обнаружение аномалий может быть весьма полезным в случае систем механической обработки. Он может точно определить моменты времени, когда станок вел себя неустойчиво, что в конечном итоге может помочь в планировании более эффективных операций обработки. Обнаружение аномалий в сочетании с прогнозированием помогает предсказать, существует ли риск отказа машины в будущем и в какой момент времени.
Просто взглянув на график повторной выборки набора данных «vib_table» для случая 11, мы увидим, что существует аномалия.
Попробуем реализовать алгоритм, способный обнаруживать эту аномальную точку.
Одним из широко используемых алгоритмов обнаружения аномалий является алгоритм изолированного леса. Это неконтролируемый алгоритм, который выявляет аномалии, изолируя выбросы в данных. Изолирующий лес основан на алгоритме дерева решений.
from sklearn.preprocessing import StandardScaler from sklearn.ensemble import IsolationForest outlier_fraction = float(0.001) scaler = StandardScaler() np_scaled = scaler.fit_transform(df5_vib_table_resampled['vib_table'].values.reshape(-1, 1)) data = pd.DataFrame(np_scaled) model = IsolationForest(contamination = outlier_fraction) model.fit(data) df5_vib_table_resampled['anomaly'] = model.predict(data) # visualization fig, ax = plt.subplots(figsize=(10,6)) a = df5_vib_table_resampled.loc[df5_vib_table_resampled['anomaly'] == -1, ['vib_table']] #anomaly ax.plot(df5_vib_table_resampled.index, df5_vib_table_resampled['vib_table'], color='black', label = 'Normal') ax.scatter(a.index,a['vib_table'], color='red', label = 'Anomaly') plt.legend() plt.show();
Немного поиграв с параметром outlier_fraction, алгоритм смог успешно обнаружить аномалию.
Рекомендации
[1] К. Гебель, Управление неопределенностью при проверке датчиков, объединение датчиков и диагностика механических систем с использованием методов мягких вычислений, доктор философии. Диссертация кафедры машиностроения Калифорнийского университета в Беркли, 1996 г.
[2] Марко Пейшейро, Прогнозирование временных рядов в Python, Manning Publishers, 2019 г.
[3] https://www.statisticshowto.com/mean-absolute-percentage-error-mape/
