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

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

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

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

Проблемы

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

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

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

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

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

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

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

Задача

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

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

Данные в виде .csv файлов, которые содержат 420 записей концентраций газов ????₂, ????????, ????₂????₄ и ????₂????₂ в трансформаторном масле, представленных в виде временной зависимости. Также доступна разметка - каждому файлу соответствует число - оставшееся время до отказа оборудования на конец файла в точках. Период между моментами времени - 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 (высокий оранжевый столбец справа на графике). Это предельное значение, выше которого остаточному ресурсу присвается константное значение 1093, то есть остаточный ресурс довольно большой, и не требуется определение точного значения с точки зрения целей диагностики.

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

Чтобы уменьшить ошибку и решить проблему моделей с непониманием и невозможностью довольно точно предсказывать значение 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, машинному обучению и искусственному интеллекту. Буду рад вашей подписке !