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

Для кого: Статья предназначена для специалистов начального/среднего уровня. Может быть полезна, например, диагностам и техническим специалистам, знакомящимся с применением машинного обучения для диагностики, дата сайентистам, желающим узнать о применении ML в промышленности, специалистам, работающим на стыке областей знаний или тем, кто только начал свой путь в машинное обучение и анализ данных!

Материалы: Код из статьи доступен на Kaggle по ссылке для самостоятельного изучения, доработки и тд. Данные с описанием доступны на Kaggle по ссылке.

1. Общая информация о задаче

Проблемы

Значительное количество силовых трансформаторов АЭС эксплуатируются с увеличенным сроком службы (назначенный срок службы - 25 лет, но не волнуйтесь - продление срока эксплуатации вполне безопасно!). Учитывая увеличенный срок службы энергоблоков АЭС, возникает необходимость контроля технического состояния силовых трансформаторов.

Текущий процесс диагностики

Для начала представим пару терминов:

Хроматография — метод разделения и анализа смесей веществ, а также изучения физико-химических свойств веществ. Основан на распределении веществ между двумя фазами — неподвижной (твёрдая фаза или жидкость, связанная на инертном носителе) и подвижной (газовая или жидкая фаза, элюент).

Хроматографический анализ - диагностический метод, который заключается в принудительном извлечении газов из масла (хроматографии), определении их качественного состава и количественного анализа.

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

Более подробно о текущем процессе диагностике можно почитать в СТО 34.01-23-003-2019 «Методические указания по техническому диагностированию развивающихся дефектов маслонаполненного высоковольтного электрооборудования по результатам анализа газов, растворенных в минеральном трансформаторном масле».

Задача

Для качественного технического обслуживания и ремонта необходимо заранее знать о возможных неисправностях, а также об остаточном ресурсе трансформаторного оборудования. Необходимо разработать модель, которая будет определять остаточный ресурс оборудования по последним 420 точкам.

Исходные данные

Данные в виде .csv файлов, которые содержат 420 записей концентраций газов H₂, CO, C₂H₄ и C₂H₂ в трансформаторном масле, представленных в виде временной зависимости. Также доступна разметка - каждому файлу соответствует число - оставшееся время до отказа оборудования на конец файла в точках. Период между моментами времени - 12 часов.

Про решение задачи остаточного ресурса

Подробно о разных подходах к решению задачи остаточного ресурса можно почитать в этой статье на Хабре. Многообразие подходов представлено на схеме из статьи:

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

2. Загрузка данных

Загрузка библиотек

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

PATH = '/kaggle/input/transformer-time/'

Загрузка разметки обучающей выборки

Мы имеем файл разметки (по сути словарь), в котором представлены данные по остаточному времени работы трансформатора до отказа: каждому файлу с концентрациями газов соответствует 1 число - оставшееся время до отказа в момент окончания файла.

y_data = pd.read_csv(PATH + 'train.csv', index_col='id')
y_data.head()

где:

• id - название файла, которые лежит в PATH с изменяемыми во времени параметрами трансформатора. 

• predicted - остаточное время работы до отказа оборудования

Загрузка обучающей выборки (параметров трансформатора)

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

X_data = {}
for row in y_data.iterrows():
    file_name = row[0]
    path = PATH + f'data_train/data_train/{file_name}'
    X_data[file_name] = pd.read_csv(path)

• Размерность каждого файла: (420, 4)

• Общее количество файлов с разметкой: 2100

Примерно так выглядят первые 2 строки каждого файла:

X_data[file_name].head(2)

На Kaggle доступна обучающая выборка с разметкой и X_test без правильных ответов (разметки). Для целей демонстрации возможности и качества решения задачи будем использовать только Обучающую выборку. Именно ее мы будем разбивать на обучающую и валидационную (отложенную) выборки для обучения и тестирования полученных моделей.

3. Моделирование

Загрузка библиотек

from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import mean_absolute_error, f1_score
from catboost import CatBoostRegressor
from tsfresh.feature_extraction import extract_features, MinimalFCParameters

Бейзлайн средним значением RUL

В качестве самого простого бейзлайн без машинного обучения возьмем среднее значение на обучающей выборке и будем им прогнозировать остаточный ресурс. Для этого реализуем следующий пайплайн:

1. Разбиение выборки на тренировочную и валидационную

2. Расчет среднего значения остаточного ресурса на обучающей выборке и построение прогноза средним

3. Оценка качества модели на валидационной выборке

В начале наглядно продемонстрируем все шаги (далее код будет под катом):

1. Для начала разделим выборку на обучающую и тестовую

y_train, y_val = train_test_split(y_data, test_size=0.25, random_state=1)
print(f'Train set shape: {y_train.shape}')
print(f'Test set shape: {y_val.shape}')

Train set shape: (1575, 1)

Test set shape: (525, 1)

2. Далее посчитаем среднее и построим прогноз

y_train_pred = [y_train.mean().values[0]] * len(y_train)
y_val_pred = [y_train.mean().values[0]] * len(y_val)

3. Рассчитаем метрики

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 219.97

MAE on test set = 218.38

Мы получили первые метрики, от которых можно отталкиваться. Если модели машинного обучения не смогут значительно улучшить метрики, то их использование не имеет смысла.

Подход на основе моделей регрессии

Давайте попробуем что-то более сложное и с машинным обучением. Реализуем следующий пайплайн для разных вариантов агрегации и моделей:

1. Агрегация данных (простая агрегация/агрегация с помощью TSFresh)

2. Разбиение выборки на тренировочную и валидационную

3. Нормализация данных (при необходимости)

4. Обучение модели (линейная регрессия/градиентный бустинг)

5. Инференс модели

6. Оценка качества модели на валидационной выборке

Подход 1: Простая агрегация (средним) + Линейная регрессия

Под простой агрегацией понимается усреднение параметров по колонкам в рамках каждого файла (4 значения средних на 4х признаках соответственно). Для каждого одномерного временного ряда получается 1 число - среднее по ряду. То есть для каждого файла мы получаем 4 признака (4 средних из 4х временных рядов).

y = y_data.copy()
X = pd.concat([X_data[file].mean() for file in y_data.index], axis=1).T
X.index = y.index
X = X.add_suffix('_mean')

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

StSc = StandardScaler()
X_train_sc = StSc.fit_transform(X_train)
X_val_sc = StSc.transform(X_val)

lr = LinearRegression()
lr.fit(X_train_sc, y_train)

y_train_pred = lr.predict(X_train_sc)
y_val_pred = lr.predict(X_val_sc)

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 177.7

MAE on test set = 175.12

Подход 2: Простая агрегация (средним) + Градиентный бустинг (CatBoost)

y = y_data.copy()
X = pd.concat([X_data[file].mean() for file in y_data.index], axis=1).T
X.index = y.index
X = X.add_suffix('_mean')

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

cbr = CatBoostRegressor(random_state=1, verbose=0)
cbr.fit(X_train, y_train)

y_train_pred = cbr.predict(X_train)
y_val_pred = cbr.predict(X_val)

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 100.6

MAE on test set = 161.67

Подход 3: Агрегация данных с помощью TSFresh + Линейная регрессия

Библиотека для агрегации и выделения признаков из временных рядов TSFresh позволяет выделять гораздо более сложные признаки по сравнению с простым средним. Мы для демонстрации будем использовать минимальные 10 признаков: sum, median, mean, length, std, variance, RMS, max, abs_max, min.

X = pd.concat([X_data[file].assign(id=file) for file in y_data.index], axis=0, ignore_index=True)
y = y_data.copy()

settings = MinimalFCParameters()
X = extract_features(X, 
                     column_id="id", 
                     default_fc_parameters=settings).loc[y.index]

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

StSc = StandardScaler()
X_train_sc = StSc.fit_transform(X_train)
X_val_sc = StSc.transform(X_val)

lr = LinearRegression()
lr.fit(X_train_sc, y_train)

y_train_pred = lr.predict(X_train_sc)
y_val_pred = lr.predict(X_val_sc)

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 140.82

MAE on test set = 148.21

Подход 4: Агрегация данных с помощью TSFresh + Градиентный бустинг (CatBoost)

X = pd.concat([X_data[file].assign(id=file) for file in y_data.index], axis=0, ignore_index=True)
y = y_data.copy()

settings = MinimalFCParameters()
X = extract_features(X, 
                     column_id="id", 
                     default_fc_parameters=settings).loc[y.index]

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

cbr = CatBoostRegressor(random_state=1, verbose=0)
cbr.fit(X_train, y_train)

y_train_pred = cbr.predict(X_train)
y_val_pred = cbr.predict(X_val)

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 37.2

MAE on test set = 91.37

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

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

plt.figure(figsize=(6, 3))
plt.hist(y_val_pred, bins=30, alpha=0.6, label='Predicted values on the test set')
plt.hist(y_val, bins=30, alpha=0.6, label='True values of the test set')
plt.legend()
plt.show()

Подход на основе ансамбля (классификация+регрессия)

Чтобы уменьшить ошибку и решить проблему моделей с непониманием и невозможностью довольно точно предсказывать значение 1093. Для этого перед прогнозом значения остаточного ресурса моделью регрессии поставим модель классификации, которая будет говорить, относится ли состояние трансформатора к одному из вариантов: "RUL=1093" или "RUL<1093". В данных нет "RUL>1093", так как данные уже предобработаны до нас. Реализуем следующий пайплайн для разных вариантов агрегации и моделей:

1. Агрегация данных с помощью TSFresh

2. Разбиение выборки на тренировочную и валидационную

3. Нормализация данных (при необходимости)

4. Обучение модели бинарной классификации (0 - "RUL<1093"; 1 - "RUL=1093")

5. Обучение модели регрессии (линейная регрессия/градиентный бустинг) на данных "RUL<1093"

6. Инференс модели бинарной классификации

7. Инференс модели регрессии на данных, где модель классификации предсказала 0

8. Оценка качества модели на валидационной выборке

Подход 1: TSFresh + Логистическая регрессия + Линейная регрессия

y = (y_data == 1093).astype(int)
X = pd.concat([X_data[file].assign(id=file) for file in y_data.index], axis=0, ignore_index=True)

settings = MinimalFCParameters()
X = extract_features(X, 
                     column_id="id", 
                     default_fc_parameters=settings).loc[y.index]

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

StSc = StandardScaler()
X_train_sc = StSc.fit_transform(X_train)
X_val_sc = StSc.transform(X_val)

LogReg = LogisticRegression()
LogReg.fit(X_train_sc, y_train)

y_train_pred_outliers = pd.DataFrame(LogReg.predict(X_train_sc), 
                                     index=y_train.index, 
                                     columns=y_train.columns)
y_val_pred_outliers = pd.DataFrame(LogReg.predict(X_val_sc), 
                                   index=y_val.index, 
                                   columns=y_val.columns)

f1_train = round(f1_score(y_train, y_train_pred_outliers), 2)
f1_val = round(f1_score(y_val, y_val_pred_outliers), 2)
print(f'F1 on train set = {f1_train}')
print(f'F1 on test set = {f1_val}')

F1 on train set = 0.67

F1 on test set = 0.64

Метрики классификации могут быть значительно улучшены!

y = y_data.copy()

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

y_train_wo_outliers = y_train[y_train != 1093].dropna()
X_train_wo_outliers = X_train.loc[y_train_wo_outliers.index]

StSc = StandardScaler()
X_train_wo_outliers_sc = StSc.fit_transform(X_train_wo_outliers)
X_train_sc = StSc.transform(X_train)
X_val_sc = StSc.transform(X_val)

lr = LinearRegression()
lr.fit(X_train_wo_outliers_sc, y_train_wo_outliers)

y_train_pred = pd.DataFrame(lr.predict(X_train_sc), 
                            index=y_train.index, 
                            columns=y_train.columns)
y_val_pred = pd.DataFrame(lr.predict(X_val_sc), 
                          index=y_val.index, 
                          columns=y_val.columns)

ind = y_train_pred_outliers[y_train_pred_outliers['predicted']==1].index
y_train_pred.loc[ind] = 1093

ind = y_val_pred_outliers[y_val_pred_outliers['predicted']==1].index
y_val_pred.loc[ind] = 1093

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 133.12

MAE on test set = 137.01

Мы, конечно, улучшили результат по сравнению с обычной моделью линейной регрессии, но этого все равно не достаточно.

Подход 2: TSFresh + Логистическая регрессия + Градиентный бустинг

Этап 1 - классификации остается прежним, код не меняется, поэтому представим код только для этапа 2:

y = y_data.copy()

X_train, X_val, y_train, y_val = train_test_split(X, y, test_size=0.25, random_state=1)

y_train_wo_outliers = y_train[y_train != 1093].dropna()
X_train_wo_outliers = X_train.loc[y_train_wo_outliers.index]

StSc = StandardScaler()
X_train_wo_outliers_sc = StSc.fit_transform(X_train_wo_outliers)
X_train_sc = StSc.transform(X_train)
X_val_sc = StSc.transform(X_val)

cbr = CatBoostRegressor(random_state=1, verbose=0)
cbr.fit(X_train_wo_outliers_sc, y_train_wo_outliers)

y_train_pred = pd.DataFrame(cbr.predict(X_train_sc), 
                            index=y_train.index, 
                            columns=y_train.columns)
y_val_pred = pd.DataFrame(cbr.predict(X_val_sc), 
                          index=y_val.index, 
                          columns=y_val.columns)

ind = y_train_pred_outliers[y_train_pred_outliers['predicted']==1].index
y_train_pred.loc[ind] = 1093

ind = y_val_pred_outliers[y_val_pred_outliers['predicted']==1].index
y_val_pred.loc[ind] = 1093

mae_train = round(mean_absolute_error(y_train, y_train_pred), 2)
mae_val = round(mean_absolute_error(y_val, y_val_pred), 2)
print(f'MAE on train set = {mae_train}')
print(f'MAE on test set = {mae_val}')

MAE on train set = 77.11

MAE on test set = 106.49

4. Заключение

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

Кстати, результат может быть улучшен в 2-3 раза (снижение ошибки) по сравнению с полученным в статье, поэтому можете пробовать! К тому же есть еще довольно много гипотез для проверки, например:

1. Расширение набора признаков средствами TSFresh (до нескольких тысяч признаков. последующим отбором)

2. Обучение моделей регрессии после классификации не на данных "RUL<1093", а на всех

3. Использование результатов модели классификации в качестве признака на вход моделям регрессии

4. Использование более сложной модели классификации для повышения качества разделения кейсов "RUL=1093" и "RUL<1093"

5. (Принципиально другой подход к постановке и решению задачи) Прогнозирование остаточного ресурса на основе прогнозирования временных рядов до пересечения с уставкой отказа (можно найти в СТО 34.01-23-003-2019)

6. и многие другие

Я создал телеграм канал DataKatser, где появляюсь гораздо чаще и делюсь своими мыслями и интересными кейсами по data science, машинному обучению и искусственному интеллекту. Буду рад вашей подписке !