Изучение набора данных НАСА по ТРДД

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

Объяснение опасностей и пример реализации

‹отказ от ответственности: я стремлюсь продемонстрировать влияние различных методов и вариантов, сделанных во время разработки модели. Эти эффекты часто демонстрируются с помощью набора тестов, что считается (очень) плохой практикой, но помогает в образовательных целях. ›

Добро пожаловать в очередной выпуск серии «Изучение набора данных НАСА по турбовентиляторным двигателям». Это будет четвертый и последний анализ первого набора данных (FD001), в котором все двигатели работают в одинаковых рабочих условиях и развивают одну и ту же неисправность.

В моем последнем посте мы углубились в анализ временных рядов и исследовали модели распределенных задержек для прогнозного обслуживания. Последняя модель показала неплохие результаты со среднеквадратичным значением 20,85. Сегодня мы рассмотрим анализ выживаемости. Техника, которую я очень хочу попробовать, поскольку я слышал и читал несколько раз, это может быть подходящим подходом для профилактического обслуживания. Однако я никогда не встречал примера реализации, которая удовлетворила бы мое любопытство. Во-первых, что такое анализ выживаемости?

Праймер для анализа выживаемости

Анализ выживаемости зародился в медицинском секторе, чтобы ответить на вопросы о продолжительности жизни конкретных групп населения. Если вы знаете чей-то возраст и можете предсказать чью-то жизнь, вы также можете оценить, сколько времени этому человеку осталось жить. Этот метод применяется, например, в эпидемиологии или исследованиях для лечения заболеваний. Однако его также можно применить во многих других случаях, когда данные состоят из продолжительности и событий, зависящих от времени, таких как прогнозирование оттока и профилактическое обслуживание. Ниже я кратко резюмирую несколько ключевых концепций, используемых в анализе выживаемости [1, 2]:

Событие: возникновение интересующего явления, в нашем случае поломка двигателя.
Продолжительность: продолжительность относится ко времени начала наблюдения. до события или остановки наблюдения
Цензура: Цензура происходит, когда наблюдения остановлены, но у интересующего субъекта еще не было своего «события».
Выживание функция: функция выживания возвращает вероятность выживания в / прошедшее время t
функция опасности: функция опасности возвращает вероятность события, произошедшего в момент времени t, при условии, что событие еще не произошло до времени t

Один из привлекательных аспектов анализа выживания для меня - это возможность включать субъектов (или в нашем случае машины) в модель, для которой еще не было своего события. В более традиционном машинном обучении вы бы исключили «неполные» или подвергнутые цензуре предметы из набора данных, что может исказить результаты [3]. Объяснив некоторые основы, пора начинать!

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

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

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

Подготовка данных

Для моделей, которые мы будем использовать позже, требуется столбец событий. Итак, давайте добавим столбец с разбивкой, показывающий, вышел ли двигатель из строя (1) или все еще работает (0).

Затем нам нужно указать время начала и окончания каждого наблюдения. Поскольку в наборе данных есть непрерывные измерения в течение временных циклов, каждое наблюдение будет составлять всего один цикл. Мы можем использовать столбец time_cycles, чтобы указать конец наблюдения, и мы добавим начальный столбец, который равен time_cycles — 1, чтобы указать начало наблюдения.

Обратите внимание на значения столбцов time_cycles, RUL, разбивки и начала, чтобы проверить, соответствует ли проведенная нами подготовка нашим ожиданиям, выглядит хорошо!

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

Кривая Каплана-Мейера

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

Этот метод уже дает нам грубый инструмент для оценки вероятности выживания в прошедшее время t для двигателя из той же совокупности. Например, у двигателей есть 100% вероятность выживания первых 128 time_cycles. После этого первые двигатели начинают выходить из строя, но все еще существует 46% вероятность того, что двигатель выживет после 200 time_cycles. Обратите внимание, что этот метод только указывает вероятность выживания после определенного момента, но не может экстраполировать за пределы предоставленных данных.

Модели пропорциональных рисков Кокса

Распространенной моделью, которая предоставляет больше информации, является модель пропорциональных опасностей Кокса. Он вычисляет коэффициент опасности, указывая, например, риск отказа, например у двигателей, работающих на определенных настройках, вероятность отказа в 1,85 раза выше, чем у двигателей, работающих на других настройках.

Реализация CoxPH пакетов жизненного цикла python также поставляется с изящным методом «pred_expectation», который дает вам прямой способ оценить время до события. Из-за этого метода pred_expectation я изо всех сил старался применить модель CoxPH к нашему набору данных. Желая использовать деградацию движка с течением времени, я использовал «cluster_col», чтобы указать движок unit_nr, пытаясь учесть в модели несколько наблюдений для каждого движка. К сожалению, результаты были довольно плохими. Преисполненный решимости добиться успеха, я обратился к автору жизненного цикла Кэмерон Дэвидсон-Пилон. Здесь я узнал, что «cluster_col» не предназначен для обозначения отсчетов, связанных со временем, а для обозначения групп с независимыми от времени наблюдениями. Например, чтобы указать разные группы обработки или группы двигателей, работающих с разными рабочими настройками.

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

CoxTimeVaryingFitter

Обучение модели довольно простое: вы создаете экземпляр модели и вызываете метод подгонки, передавая набор данных, id_col для указания уникальных механизмов, event_col, чтобы указать, произошла ли поломка, и столбцы start и stop, чтобы модель могла интерпретировать продолжительность наблюдения.

# returns
# Iteration 8: norm_delta = 0.00000, step_size = 1.00000, ll = -63.86010, newton_decrement = 0.00000, seconds_since_start = 0.1
# Convergence completed after 8 iterations.
# <lifelines.CoxTimeVaryingFitter: fitted with 18627 periods, 100 subjects, 54 events>

Далее осмотрим модель.

# returns
# <lifelines.CoxTimeVaryingFitter: fitted with 18627 periods, 100 subjects, 54 events>
#          event col = 'breakdown'
# number of subjects = 100
#  number of periods = 18627
#   number of events = 54
#     log-likelihood = -63.86
#   time fit was run = 2020-09-02 15:26:11 UTC

Глядя на сводку модели, нас интересуют логарифмическая вероятность, p-значения и exp (coef). Логарифмическое правдоподобие показывает степень соответствия, но только по сравнению с другими аналогичными моделями, имеющими меньшее количество функций. Логарифмическое правдоподобие, близкое к 0, считается лучшим (не путать с логарифмическим отношением правдоподобия!). Если посмотреть на значения p, то значения для датчиков 9 и 15 довольно велики при p ›0,50. Однако удаление датчиков 9 и 15 вернуло логарифмическую правдоподобность -64,20, таким образом не улучшив точность подгонки [4, 5].

Эксп (коэф) показывает риск опасности масштабирования. Например. приращение на 1 единицу значений датчика датчика 11 увеличивает риск поломки на 167,43 [6]. График по существу отображает коэффициенты и доверительные интервалы признаков.

Прогнозировать и оценить

Когда модель обучена, пора приступить к оценке. Риск отказа (или опасности) зависит от базовой и частичной опасности (см. Формулу ниже).

Поскольку наши двигатели относятся к однородной группе (например, все двигатели работают в одинаковых рабочих условиях), их базовая опасность одинакова. Поэтому нам нужно только проверить частичную или частичную логарифмическую опасность, чтобы получить указание на риск отказа. Частичная опасность имеет значение только по отношению к другим частичным опасностям, исходящим от той же группы населения. Вероятность поломки двигателя с частичной опасностью 2e⁶ в два раза выше, чем у двигателя с частичной опасностью 1e⁶. Поскольку значения частичных опасностей довольно велики, легче отобразить журнал частичных опасностей. Однако частичная опасность журнала снижает интерпретируемость. На каждые 1 единицу увеличения логарифмической частичной опасности одного двигателя по сравнению с другим вероятность поломки становится в 2,718 (или e) раз больше.

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

Сравнивая log_partial_hazard с вычисленным RUL, вы можете увидеть, что в целом он довольно хорошо сообщает о неизбежности поломки (здесь показаны первые 10). Более высокие значения log_partial_hazards возвращаются для двигателей, более подверженных риску выхода из строя. Однако это не всегда точно, например, опасность для двигателя 16 намного выше, чем опасность для двигателя 15, хотя двигатель 15 выходит из строя раньше.

Построение графика всех log_partial_hazards против вычисленного RUL дает следующий график с четким видимым трендом.

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

Примечание: здесь можно установить пороговое значение для log_partial_hazard, после которого следует проводить техническое обслуживание. Однако, как обсуждалось ранее, это на самом деле не информирует вас о ПРАВИЛЕ. Вы можете разработать модель временных рядов, чтобы предсказать, когда будет достигнут этот порог, чтобы получить более точный прогноз «времени до события».

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

Регрессия частичной логарифмической опасности в RUL

Сначала мы спрогнозируем log_partial_hazard для каждого наблюдения в цензурированном обучающем наборе и проверим его график разброса.

Вы можете ясно видеть влияние нашего клиппера RUL в верхней части графика, но разброс был бы еще больше без отсечения. Взаимосвязь, отображаемая на диаграмме рассеяния, является нелинейной, возможно, экспоненциальной. Разброс довольно небольшой, что затрудняет определение RUL на основе log_partial_hazard, но давайте посмотрим, как мы будем жить.

Мы подберем экспоненциальную модель, чтобы вывести RUL из log_partial_hazard, используя scipy curve_fit [8]. После вывода RUL мы оценим его по сравнению с вычисленным RUL для обучающего и тестового набора, чтобы получить представление о его точности.

Сначала мы определяем функцию оценки.

Затем определяется и аппроксимируется экспоненциальная модель с использованием кривой Сципи.

# returns
# array([8.85954699e+01, 4.35302167e-02])

Наконец, готовится набор тестов, и оцениваются предсказания поездов и тестов.

# returns
# train set RMSE:26.302858422243872, R2:0.5487597218187095
# test set RMSE:27.135244169256758, R2:0.5736091039386049

Среднеквадратичное значение 27,13 уже на 15% больше по сравнению с нашей базовой моделью, которая имела среднеквадратичное значение 31,95. Теперь давайте потренируемся на полном наборе данных и посмотрим, как работает модель.

Повторить для полного набора данных

# returns
# Iteration 8: norm_delta = 0.00000, step_size = 1.00000, ll = -114.77106, newton_decrement = 0.00000, seconds_since_start = 0.2 
# Convergence completed after 8 iterations.
# <lifelines.CoxTimeVaryingFitter: fitted with 20631 periods, 100 subjects, 100 events>
# train set RMSE:26.226364780597272, R2:0.6039289060308352
# test set RMSE:26.580988808221285, R2:0.5908498441209542

Окончательное среднеквадратичное значение составляет 26,58, что на 16,8% больше нашего базового уровня, но не приближается к решениям для УВО (RMSE = 20,54) или анализа временных рядов (RMSE = 20,85). Частично неточность может быть объяснена путем подбора другой модели поверх предсказанного log_partial_hazard, что приводит к ошибкам поверх ошибок (поскольку ни одна модель не является идеальной).

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

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

В целом, я считаю эту технику достаточно интересной, и не помешало бы узнать о ней немного больше!

Я хотел бы выразить особую благодарность автору жизненных линий Кэмерон Дэвидсон-Пилон за то, что он нашел время и дал мне несколько советов о том, как наилучшим образом использовать пакет жизненных линий для имеющегося набора данных. Кроме того, я хотел бы поблагодарить Wisse Smit и Maikel Grobbe за их вклад и рецензирование моей статьи. Вы можете найти полную записную книжку на моей странице на github здесь.

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

Ссылки:
[1] https://lifelines.readthedocs.io/en/latest/Survival%20Analysis%20intro.html
[2] https://en.wikipedia.org/ wiki / Survival_analysis
[3] https://lifelines.readthedocs.io/en/latest/Survival%20Analysis%20intro.html#censoring
[4] https: // stats. idre.ucla.edu/other/mult-pkg/faq/general/faqhow-are-the-likelihood-ratio-wald-and-lagrange-multiplier-score-tests-different-andor-similar/
[ 5] https://www.reddit.com/r/statistics/comments/23sk6h/what_does_a_loglikelihood_value_indicate_and_how/
[6] https://medium.com/@zachary.james.angell/applying-survival -analysis-to-customer-churn-40b5a809b05a
[7] https://lifelines.readthedocs.io/en/latest/Time%20varying%20survival%20regression.html
[8] Https://stackoverflow.com/questions/52930401/how-to-get-a-robust-nonlinear-regression-fit-using-scipy-optimize-least-squares