Создание seq2seq модели для задачи регрессии¶
В этом примере мы обучим глубокую нейросеть предсказывать остаточный ресурс газотурбинного двигателя.
Загрузка и распаковка данных¶
В открытом датасете CMAPSSData хранятся имитационные данные о работе двигателей: 100 учебных примеров и 100 примеров для валидации. Обучающие данные содержат временные ряды, состоящие из последовательности измеренных параметров, записанных от пуска до наступления отказа. Данные для валидации обрываются до наступления события отказа.
В каждой строке содержится 26 переменных:
- 1: Номер двигателя
- 2: Наработка (в циклах)
- 3–5: Настройки двигателя
- 6–26: Измерения датчиков 1–21
Подготовим рабочее окружение и загрузим набор данных.
using DelimitedFiles
dataTrain = Float32.( readdlm( "CMAPSSData/train_FD001.txt" ));
Подготовим обучающие данные¶
В результате этой обработки мы получим массивы XTrain и YTrain, содержащих временные ряды, состоящие из признаков (предикторов, predictors) и целевых переменных (таргетов, targets).
function dataPreparation( dataTable )
numObservations = maximum( dataTrain[:,1] )
predictors = []
responses = []
for i ∈ 1:numObservations
idx = dataTable[:,1] .== i
push!( predictors, dataTable[idx, 3:end] )
timeSteps = dataTable[idx,2]
push!( responses, reverse(timeSteps) )
end
return predictors, responses
end
XTrain, YTrain = dataPreparation( dataTrain );
Удалить константные признаки¶
Неизменные признаки ничего не добавляют нашему процессу обучения, поэтому для ускорения вычислений мы от них избавимся. А именно – удалим из датасета столбцы, у которых минимум и максимум являются одинаковыми значениями.
m = [ minimum( table, dims=1 ) for table in XTrain ];
M = [ maximum( table, dims=1 ) for table in XTrain ];
idxConstant = [ m[i] .== M[i] for i in 1:length(m) ];
constant_features = vec( maximum( vcat(idxConstant...), dims=1 ) );
XTrain = [ xtrain[:, .!constant_features] for xtrain in XTrain ];
numFeatures = size( XTrain[1], 2 )
Осталось 16 признаков для каждого наблюдения (стоит помнить, что все наблюдения имеют различную длину, каждый двигатель наработал разное количество циклов до отказа).
Нормализация значений признаков¶
Еще один прием, который обычно ускоряет обучение, заключается в том, чтобы подготовить признаки, задав каждой серии математическое ожидание 0 и дисперсию 1. Соберём сигналы со всех двигателей в один столбец.
using Statistics
mu = vec(mean( vcat([xtrain for xtrain in XTrain]...), dims=1 ));
sig = vec(std( vcat([xtrain for xtrain in XTrain]...), dims=1 ));
XTrain = [ mapslices(row -> (row .- mu) ./ sig, xtrain; dims=2) for xtrain in XTrain ];
gr()
plot(
heatmap( XTrain[1], cbar=false, yflip=true, title="Наблюдения (predictors)", ylabel="Время (циклы)" ),
heatmap( reshape(YTrain[1],:,1), yflip=true, title="Целевая переменная: остаток ресурса (target)" ),
size=(800,400)
)
plot!( titlefont = font(9) )
Ограничение для выходной переменной¶
Мы будем предсказывать одно число – количество циклов работы, остающихся до наступления отказа. Нам важнее, чтобы предсказания были более точными ближе к концу жизненного цикла.
Поэтому ограничим сверху предсказываемую переменную ресурсом в 150 циклов наработки на отказ.
thr = 150
YTrain = [ min.(ytrain, 150) for ytrain in YTrain ];
Так выглядят показания датчиков для некоторого двигателя в датасете.
obs_id = 77
dl = size( XTrain[obs_id], 1 );
cl = maximum( findall( YTrain[obs_id] .== thr ))
a = plot( XTrain[obs_id], leg=:false, title="Показания датчиков", label=:none )
vline!(a, [dl], lw=3, lc=:red, label=:none )
plot!(a, [1, cl], [-4,-4], fillrange=[4,4], fillcolor=:springgreen, fillalpha=0.3, linealpha=0.0, label=:none )
b = plot( YTrain[obs_id], title="Наработка на отказ (в циклах)", label="Ресурс" )
vline!(b, [dl], lw=3, lc=:red, label="Отказ" )
plot!(b, [1, cl], [0,0], fillrange=[thr+20,thr+20], fillcolor=:springgreen, fillalpha=0.3, linealpha=0.0, label="Полный ресурс" )
plot(a, b, layout=(2,1) )
Отсортируем данные¶
Это необязательный этап, когда-то он поможет нам разбить данные на более оптимальные батчи, но для начала нам просто нужно узнать, чему равна длина наибольшей последовательности в выборке.
# Упорядочить данные по размеру
sequenceLength = [ size(xtrain,1) for xtrain in XTrain ];
idx = sortperm( sequenceLength, rev=true );
XTrain = XTrain[idx];
YTrain = YTrain[idx];
bar( sort(sequenceLength), xlabel="Последовательность", ylabel="Длительность", title="Отсортированные данные", xflip=true)
Архитектура нейросети¶
Мы обучим две рекуррентные нейросети и сравним результаты:
- RNN нейросеть (обычный рекуррентный слой, затем несколько полносвязанных)
- LSTM нейросеть (вместо базового рекуррентного слоя чуть более сложный, LSTM слой, остальные слои те же самые)
До обновления библиотеки
Fluxдо последней версии следующий код будет выдавать безвредное сообщение об ошибкеError during loading of extension DiffEqBaseZygoteExt..., которое не мешает работе. Просто нажмите на него правой кнопкой и выберитеУдалить выбранное.
Pkg.add( "Flux" ); # Установка библиотеки на случай ее отсутствия
using Flux
numResponses = size( YTrain[1], 2 )
numHiddenUnits = 100
numHiddenUnits2 = 40
rnn = Chain( RNN(numFeatures => numHiddenUnits),
Dense(numHiddenUnits => numHiddenUnits2),
Dropout( 0.5 ),
Dense(numHiddenUnits2 => numResponses) ) |> f32;
lstm = Chain(
LSTM(numFeatures => numHiddenUnits),
Dense( numHiddenUnits => 50 ),
Dropout( 0.5 ),
Dense( 50 => numResponses )
) |> f32;
Запомним эти две нейросети в том состоянии, которое они имели до начала обучения, чтобы потом сравнить их с обученными сетями.
rnn_empty = deepcopy( rnn );
lstm_empty = deepcopy( lstm );
Функция потерь (RMSE) сформулирована следующим образом:
function loss(x, y)
s = 0
for (xi,yi) in zip(x,y)
Flux.reset!( model );
xi_pred = [ model(xi[i,:])[end] for i in 1:size(xi,1) ];
s = s + (mse(xi_pred, yi)) .* 2
end
return sqrt( s )
end
Она возвращает нам сумму ошибок всех прогнозов для каждой последовательности в выборке.
Циклы обучения¶
Сперва обучим RNN нейросеть.
epochs = 100;
Время обучения этой сети до 100-й эпохи составит 5-7 минут.
model = rnn
opt = Flux.Optimiser(ClipValue(1), ADAM(0.01))
θ = Flux.params( model )
for epoch ∈ 1:epochs
print( epoch, " ")
∇ = gradient(θ) do
loss( XTrain, YTrain )
end
Flux.update!(opt, θ, ∇)
end
# Сохраним полученную нейросеть в файл
using JLD2
if !isfile( "rnn.jld2" ) jldsave("rnn.jld2"; rnn); end
Для сравнения обучим LSTM нейросеть (осторожно, время обучения до 100-й эпохи превышает 15 минут). Добавим также что количество нейронов в скрытых слоях, конечно, влияет на время обучения, но отнюдь не так сильно, как количество эпох.
model = lstm
opt = Flux.Optimiser(ClipValue(1), ADAM(0.01))
θ = Flux.params( model )
for epoch ∈ 1:epochs
print( epoch, " ")
∇ = gradient(θ) do
loss( XTrain, YTrain )
end
Flux.update!(opt, θ, ∇)
end
# Сохраняем нейросеть
if !isfile( "lstm.jld2" ) jldsave("lstm.jld2"; lstm); end
LSTM нейросети могут научиться учитывать более долгосрочные зависимости, чем RNN. В обычных рекуррентных нейросетях влияние последних измерений на результат всегда выше, чем влияние более ранних измерений, то есть RNN имеют тенденцию "забывать далекое прошлое".
Но, из-за большего количества уравнений в каждом нейроне, LSTM нейросеть училась в 2-3 раза медленнее, чем RNN.
Проверка результатов¶
Нейросети довольно хорошо предсказывают тренировочные данные:
test_id = 6
xi = XTrain[test_id]
yi = YTrain[test_id]
# Прогнозы до обучения
y_pred_rnn_empty = [ rnn_empty(xi[i,:])[end] for i in 1:size(xi,1) ];
y_pred_lstm_empty = [ lstm_empty(xi[i,:])[end] for i in 1:size(xi,1) ];
# Прогнозы после обучения
y_pred_rnn = [ rnn(xi[i,:])[end] for i in 1:size(xi,1) ];
y_pred_lstm = [ lstm(xi[i,:])[end] for i in 1:size(xi,1) ];
plot(
plot( [y_pred_rnn_empty y_pred_lstm_empty], label=["rnn" "lstm"], title="Прогнозы до обучения", lw=2 ),
plot( [y_pred_rnn y_pred_lstm yi], label=["rnn" "lstm" "истина"], title="Прогнозы обученной нейросети", lw=2 ),
size=(800,300), titlefont=font(11)
)
Как мы убедились, до обучения нейросеть всегда выдавала прогнозы, приближенные к нулю. После обучения она научилась связывать вектор из 16 измерений с скалярным прогнозом остаточного ресурса. Можно предположить, что модель учитывает предысторию процесса, поскольку в прогнозе всегда есть тренд на уменьшение, несмотря на довольно хаотичное поведение сигнала.
Изучим прогнозы по тренировочной выборке.
function get_plot(id)
xi = XTrain[id]
yi = YTrain[id]
Flux.reset!( rnn ) # Сбросим состояние нейросети, чтобы накопленная информация
Flux.reset!( lstm ) # от прошлых прогонов не повлияла на будущие прогнозы
y_pred_rnn = [ rnn(xi[i,:])[end] for i in 1:size(xi,1) ];
y_pred_lstm = [ lstm(xi[i,:])[end] for i in 1:size(xi,1) ];
return plot( [y_pred_rnn y_pred_lstm yi], legend=false, lw=2, title="Прогнозы для выборки #$id" )
end
plot( get_plot.(1:16)..., titlefont=font(7), size=(800,600) )
Можно заметить, что на тестовых данных LSTM нейросеть научилась хорошо определять начало ресурса и отслеживать динамику до конца. Но нелинейность выученной характеристики чаще проявляется в результатах обучения RNN нейросети (горб посередине функции).
Подготовим тестовые данные:
dataTest = Float32.( readdlm( "CMAPSSData/test_FD001.txt" ));
XTest, YTest = dataPreparation( dataTest );
XTest = [ xtest[:, .!constant_features] for xtest in XTest ];
mu = vec(mean( vcat([xtest for xtest in XTest]...), dims=1 ));
sig = vec(std( vcat([xtest for xtest in XTest]...), dims=1 ));
XTest = [ mapslices(row -> (row .- mu) ./ sig, xtest; dims=2) for xtest in XTest ];
# Дополнение: мы знаем остаточный ресурс двигателей в тестовой выборке
YTest = Float32.( readdlm( "CMAPSSData/RUL_FD001.txt" ));
YTest = [ collect((size(XTest[i],1)+(YTest[i])-1):-1:YTest[i]) for i in 1:length(XTest) ];
YTest = [ min.(YTest[i], 150) for i in 1:size(YTest,1) ];
Посмотрим, как модели предсказывают наработку на отказ для одной серии измерений из общей выборки (на тестовой выборке).
function get_plot( id )
xi = XTest[id]
yi = YTest[id]
Flux.reset!( rnn ) # Сбросим состояние нейросети, чтобы накопленная информация
Flux.reset!( lstm ) # от прошлых прогонов не повлияла на будущие прогнозы
y_pred_rnn = [ rnn(xi[i,:])[end] for i in 1:size(xi,1) ];
y_pred_lstm = [ lstm(xi[i,:])[end] for i in 1:size(xi,1) ];
return plot( [y_pred_rnn y_pred_lstm yi], legend=false, lw=2, title="Прогнозы для выборки #$id", c=[4 6 3] )
end
plot( get_plot.(1:16)..., titlefont=font(7), size=(800,600), legendfont=font(8) )
На тестовых данных нейросети работают чуть хуже, но именно эта информация позволяет нам выбрать дальнейшее направление для доработки архитектуры и процесса обучения.
model = rnn; println( "Качество прогноза RNN нейросети: ", loss( XTrain, YTrain ))
model = lstm; println( "Качество прогноза LSTM нейросети: ", loss( XTrain, YTrain ))
model = rnn; println( "Качество прогноза RNN нейросети: ", loss( XTest, YTest ))
model = lstm; println( "Качество прогноза LSTM нейросети: ", loss( XTest, YTest ))
Как и следовало ожидать, качество прогноза обеих нейросетей упало при переходе на тестовые данные (в 2.1 раз для LSTM и 2.7 раз для RNN).
Напоследок, изучим гистограмму ошибок прогноза ресурса на конец каждого из изучаемых сегментов жизненного цикла:
dataset_final_rul = [ ytest[end] for ytest in YTest ]
rnn_final_prediction = [ [(Flux.reset!( rnn ); rnn(xi[i,:])[end]) for i in 1:size(xi,1)][end] for xi in XTest]
lstm_final_prediction = [ [(Flux.reset!( lstm ); lstm(xi[i,:])[end]) for i in 1:size(xi,1)][end] for xi in XTest]
plot(
histogram( dataset_final_rul .- rnn_final_prediction, leg=:none, title="Ошибка прогноза RNN" ),
histogram( dataset_final_rul .- lstm_final_prediction, leg=:none, title="Ошибка прогноза LSTM" ),
size=(900,300), titlefont=font(9)
)
Можно утверждать, что средняя ошибка LSTM сети ближе к 0, чем у RNN.
Мы сохранили нейросети в файлы для дальнейшего использования.
Для их загрузки в память нужно будет сперва загрузить библиотеку
Flux.
Заключение¶
Мы обучили нейросеть предсказывать наработку на отказ и теперь можем собрать датчик для косвенных измерений чтобы отработать его в модельном окружении Engee, а также сгенерировать код для встраиваемой платформы и проводить испытания в реальных условиях.
Заметим также, что в учебных целях, ради ускорения процесса обучения, мы реализовали совсем небольшую нейросеть и обучали ее в течение довольно небольшого количества циклов. Обучение обеих моделей заняло примерно по 10 минут в 4 потока на двух процессорах (без GPU).