Загрузите набор данных

Этот пример использует Турбовентиляторный набор данных симуляции ухудшения Engine (C-MAPSS) [2]. Zip-файл содержит данные timeseries запуска к отказу для четырех различных наборов (а именно, FD001, FD002, FD003, FD004) симулированный под различными комбинациями операционных условий и режимов отказа.

Этот пример использует только набор данных FD001, который далее разделен на тестовые подмножества и обучение. Учебное подмножество содержит симулированные данные временных рядов для 100 механизмов. Каждый механизм имеет несколько датчиков, значения которых зарегистрированы в приведенном примере в непрерывном процессе. Следовательно последовательность записанных данных варьируется по длине и соответствует полному экземпляру запуска к отказу (RTF). Тестовое подмножество содержит 100 частичных последовательностей и соответствующие значения остающегося срока полезного использования в конце каждой последовательности.

Загрузите набор данных Turbofan Engine Degradation Simulation на файл с именем “CMAPSSData.zip” и разархивируйте его к папке под названием “data” в текущем каталоге.

filename = "CMAPSSData.zip";
if ~exist(filename,'file')
    url = "https://ti.arc.nasa.gov/c/6/";
    websave(filename,url);
end

dataFolder = "data";
if ~exist(dataFolder,'dir')
    mkdir(dataFolder);
end
unzip(filename,dataFolder)

Папка данных теперь содержит текстовые файлы с 26 столбцами чисел, разделенных пробелами. Каждая строка является снимком состояния данных, взятых во время одного рабочего цикла, и каждый столбец представляет различную переменную:

Предварительно обработайте обучающие данные

Загрузите данные с помощью функционального localLoadData. Функция извлекает данные из filenamePredictors и возвращает таблицу, которая содержит учебные предикторы и соответствующий ответ (i.e., RUL) последовательности. Каждая строка представляет различный механизм.

filenameTrainPredictors = fullfile(dataFolder,"train_FD001.txt");
rawTrain = localLoadData(filenameTrainPredictors);

Исследуйте данные запуска к отказу на один из механизмов.

head(rawTrain.X{1},8)

ans=8×26 table
    id    timeStamp    op_setting_1    op_setting_2    op_setting_3    sensor_1    sensor_2    sensor_3    sensor_4    sensor_5    sensor_6    sensor_7    sensor_8    sensor_9    sensor_10    sensor_11    sensor_12    sensor_13    sensor_14    sensor_15    sensor_16    sensor_17    sensor_18    sensor_19    sensor_20    sensor_21
    __    _________    ____________    ____________    ____________    ________    ________    ________    ________    ________    ________    ________    ________    ________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________

    1         1          -0.0007         -0.0004           100          518.67      641.82      1589.7      1400.6      14.62       21.61       554.36      2388.1      9046.2        1.3         47.47       521.66         2388       8138.6       8.4195        0.03          392         2388          100         39.06       23.419  
    1         2           0.0019         -0.0003           100          518.67      642.15      1591.8      1403.1      14.62       21.61       553.75        2388      9044.1        1.3         47.49       522.28       2388.1       8131.5       8.4318        0.03          392         2388          100            39       23.424  
    1         3          -0.0043          0.0003           100          518.67      642.35        1588      1404.2      14.62       21.61       554.26      2388.1      9052.9        1.3         47.27       522.42         2388       8133.2       8.4178        0.03          390         2388          100         38.95       23.344  
    1         4           0.0007               0           100          518.67      642.35      1582.8      1401.9      14.62       21.61       554.45      2388.1      9049.5        1.3         47.13       522.86       2388.1       8133.8       8.3682        0.03          392         2388          100         38.88       23.374  
    1         5          -0.0019         -0.0002           100          518.67      642.37      1582.8      1406.2      14.62       21.61          554      2388.1      9055.1        1.3         47.28       522.19         2388       8133.8       8.4294        0.03          393         2388          100          38.9       23.404  
    1         6          -0.0043         -0.0001           100          518.67       642.1      1584.5      1398.4      14.62       21.61       554.67        2388      9049.7        1.3         47.16       521.68         2388       8132.9       8.4108        0.03          391         2388          100         38.98       23.367  
    1         7            0.001          0.0001           100          518.67      642.48      1592.3      1397.8      14.62       21.61       554.34        2388      9059.1        1.3         47.36       522.32         2388       8132.3       8.3974        0.03          392         2388          100          39.1       23.377  
    1         8          -0.0034          0.0003           100          518.67      642.56        1583        1401      14.62       21.61       553.85        2388      9040.8        1.3         47.24       522.47         2388       8131.1       8.4076        0.03          391         2388          100         38.97       23.311  

Исследуйте данные об ответе на один из механизмов.

rawTrain.Y{1}(1:8)

ans = 8×1

   191
   190
   189
   188
   187
   186
   185
   184

Визуализируйте данные timeseries для некоторых предикторов.

stackedplot(rawTrain.X{1},[3,5,6,7,8,15,16,24],XVariable='timeStamp')

Удалите функции с меньшим количеством изменчивости

Функции, которые остаются постоянными навсегда, продвигаются, может негативно повлиять на обучение. Используйте prognosability функционируйте, чтобы измерить изменчивость функций при отказе.

prog = prognosability(rawTrain.X,"timeStamp");

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

idxToRemove = prog.Variables==0 | isnan(prog.Variables);
featToRetain = prog.Properties.VariableNames(~idxToRemove);
for i = 1:height(rawTrain)
    rawTrain.X{i} = rawTrain.X{i}{:,featToRetain};
end

Нормируйте учебные предикторы

Нормируйте учебные предикторы, чтобы иметь нулевое среднее значение и модульное отклонение.

[~,Xmu,Xsigma] = zscore(vertcat(rawTrain.X{:}));
preTrain = table();
for i = 1:numel(rawTrain.X)
    preTrain.X{i} = (rawTrain.X{i} - Xmu) ./ Xsigma;
end

Отсеките ответы

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

clipResponses = true;

if clipResponses
    rulThreshold = 150;
    for i = 1:numel(rawTrain.Y)
        preTrain.Y{i} = min(rawTrain.Y{i},rulThreshold);
    end
end

Этот рисунок показывает первое наблюдение и соответствующий отсеченный ответ.

Подготовка данных для дополнения

Чтобы минимизировать объем дополнения добавленного к мини-пакетам, отсортируйте обучающие данные по длине последовательности. Затем выберите мини-пакетный размер, который делит обучающие данные равномерно и уменьшает объем дополнения в мини-пакетах.

Сортировка обучающих данных длиной последовательности.

for i = 1:size(preTrain,1)
    preTrain.X{i} = preTrain.X{i}';    %Transpose training data to have features in the first dimension
    preTrain.Y{i} = preTrain.Y{i}';    %Transpose responses corresponding to the training data
    sequence = preTrain.X{i};
    sequenceLengths(i) = size(sequence,2); 
end

[sequenceLengths,idx] = sort(sequenceLengths,'descend');
XTrain = preTrain.X(idx);
YTrain = preTrain.Y(idx);

Сетевая архитектура

Глубокая архитектура сверточной нейронной сети, используемая для оценки RUL, описана в [1].

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

Выбранная сетевая архитектура применяется 1D свертка вдоль направления последовательности времени только. Это подразумевает, что порядок функций не повлияет на учебное, и только тренды в одной функции за один раз рассматриваются.

Определить сетевую архитектуру. Создайте CNN, который состоит из пяти последовательных наборов свертки 1d, нормализация партии. и, relu слой, с увеличением filterSize и numFilters как первые два входных параметра к convolution1dLayer, сопровождаемый полносвязным слоем размера numHiddenUnits и слой уволенного с вероятностью уволенного 0,5. Поскольку сеть предсказывает остающийся срок полезного использования (RUL) турбовентиляторного механизма, установите numResponses к 1 во втором полносвязном слое и слое регрессии как последний слой сети.

Чтобы компенсировать различные последовательности времени в обучающих данных, используйте Padding="causal" как входной параметр Пары "имя-значение" в convolution1dLayer.

numFeatures = size(XTrain{1},1);
numHiddenUnits = 100;
numResponses = 1;

layers = [
    sequenceInputLayer(numFeatures)
    convolution1dLayer(5,32,Padding="causal")
    batchNormalizationLayer()
    reluLayer()
    convolution1dLayer(7,64,Padding="causal")
    batchNormalizationLayer
    reluLayer()
    convolution1dLayer(11,128,Padding="causal")
    batchNormalizationLayer
    reluLayer()
    convolution1dLayer(13,256,Padding="causal")
    batchNormalizationLayer
    reluLayer()
    convolution1dLayer(15,512,Padding="causal")
    batchNormalizationLayer
    reluLayer()
    fullyConnectedLayer(numHiddenUnits)
    reluLayer()
    dropoutLayer(0.5)
    fullyConnectedLayer(numResponses)
    regressionLayer()];

Обучите сеть

Задайте trainingOptions (Deep Learning Toolbox). Обучите в течение 30 эпох с мини-пакетами размера 20 использований ‘adam’ решателя. Установите LearnRateSchedule to piecewise. Задайте скорость обучения 0.01. Чтобы препятствовать тому, чтобы градиенты взорвались, установите порог градиента к 1. Чтобы сохранить последовательности отсортированными по длине, установите 'Shuffle' к 'never'. Включите график процесса обучения и выключите командное окно выход (Verbose).

maxEpochs = 30;
miniBatchSize = 20;

options = trainingOptions('adam',...
    LearnRateSchedule='piecewise',...
    MaxEpochs=maxEpochs,...
    MiniBatchSize=miniBatchSize,...
    InitialLearnRate=0.01,...
    GradientThreshold=1,...
    Shuffle='never',...
    Plots='training-progress',...
    Verbose=0);

Обучите сеть с помощью trainNetwork.

net = trainNetwork(XTrain,YTrain,layers,options);

Постройте график слоев сети, чтобы визуализировать базовую сетевую архитектуру.

figure;
lgraph = layerGraph(net.Layers);
plot(lgraph)

Протестируйте сеть

Тестовые данные содержат 100 частичных последовательностей и соответствующие значения остающегося срока полезного использования в конце каждой последовательности.

filenameTestPredictors = fullfile(dataFolder,'test_FD001.txt');
filenameTestResponses = fullfile(dataFolder,'RUL_FD001.txt');
dataTest = localLoadData(filenameTestPredictors,filenameTestResponses);

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

for i = 1:numel(dataTest.X)
    dataTest.X{i} = dataTest.X{i}{:,featToRetain};
    dataTest.X{i} = (dataTest.X{i} - Xmu) ./ Xsigma;
    if clipResponses
        dataTest.Y{i} = min(dataTest.Y{i},rulThreshold);
    end
end

Составьте таблицу для хранения предсказанного ответа (YPred) наряду с истинным ответом (Y). Сделайте предсказания на тестовых данных с помощью predict. Чтобы препятствовать тому, чтобы функция добавила дополнение в тестовые данные, задайте мини-пакетный размер 1.

predictions = table(Size=[height(dataTest) 2],VariableTypes=["cell","cell"],VariableNames=["Y","YPred"]);

for i=1:height(dataTest)
    unit = dataTest.X{i}';   
    predictions.Y{i} = dataTest.Y{i}';
    predictions.YPred{i} = predict(net,unit,MiniBatchSize=1);
end

Показатели производительности

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

for i = 1:size(predictions,1)
    predictions.RMSE(i) = sqrt(mean((predictions.Y{i} - predictions.YPred{i}).^2));
end

Следующая гистограмма помогает в визуализации распределения значений RMSE через все тестовые механизмы.

figure;
histogram(predictions.RMSE,NumBins=10);
title("RMSE ( Mean: " + round(mean(predictions.RMSE),2) + " , StDev: " + round(std(predictions.RMSE),2) + " )");
ylabel('Frequency');
xlabel('RMSE');

Кроме того, чтобы видеть, как сетевой предиктор выполняет в данной последовательности данных в тестовых механизмах. Используйте localLambdaPlot функционируйте, чтобы построить предсказанный RUL против истинного RUL любого тестового механизма.

figure;
localLambdaPlot(predictions,"random");

Результат показывает, что архитектура глубокого обучения CNN для оценки RUL турбо данных о механизме является альтернативным подходом, чтобы предсказать RUL. Значения RMSE во всех метках времени указывают, что сеть смогла выполнить хорошо к концу данных тестовых данных о последовательности. Это предполагает, что важно иметь некоторую маленькую историю значений датчика при попытке предсказать RUL.

Функции помощника

Load Data function

Эта функция загружает данные запуска к отказу из предоставленного текстового файла, данные timeseries групп и его соответствующие значения RUL в таблице как предикторы и ответы.

function data = localLoadData(filenamePredictors,varargin)

if isempty(varargin)
    filenameResponses = []; 
else
    filenameResponses = varargin{:};
end

%% Load the text file as a table
rawData = readtable(filenamePredictors);

% Add variable names to the table
VarNames = {...
    'id', 'timeStamp', 'op_setting_1', 'op_setting_2', 'op_setting_3', ...
    'sensor_1', 'sensor_2', 'sensor_3', 'sensor_4', 'sensor_5', ...
    'sensor_6', 'sensor_7', 'sensor_8', 'sensor_9', 'sensor_10', ...
    'sensor_11', 'sensor_12', 'sensor_13', 'sensor_14', 'sensor_15', ...
    'sensor_16', 'sensor_17', 'sensor_18', 'sensor_19', 'sensor_20', ...
    'sensor_21'};
rawData.Properties.VariableNames = VarNames;

if ~isempty(filenameResponses)
    RULTest = dlmread(filenameResponses);
end

% Split the signals for each unit ID
IDs = rawData{:,1};
nID = unique(IDs);
numObservations = numel(nID);

% initialize a table for storing data
data = table(Size=[numObservations 2],...
    VariableTypes={'cell','cell'},...
    VariableNames={'X','Y'});

for i=1:numObservations
    idx = IDs == nID(i);
    data.X{i} = rawData(idx,:);
    if isempty(filenameResponses)
        % calculate RUL from time column for train data
        data.Y{i} = flipud(rawData.timeStamp(idx))-1;
    else
        % use RUL values from filenameResponses for test data
        sequenceLength = sum(idx);
        endRUL = RULTest(i);
        data.Y{i} = [endRUL+sequenceLength-1:-1:endRUL]'; %#ok<NBRAK> 
    end
end
end

Lambda Plot function

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

function localLambdaPlot(predictions,lambdaCase)

if isnumeric(lambdaCase)
    idx = lambdaCase;
else
    switch lambdaCase
        case {"Random","random","r"}
            idx = randperm(height(predictions),1); %Randomly choose a test case to plot
        case {"Best","best","b"}
            idx = find(predictions.RMSE == min(predictions.RMSE)); %Best case
        case {"Worst","worst","w"}
            idx = find(predictions.RMSE == max(predictions.RMSE)); %Worst case
        case {"Average","average","a"}
            err = abs(predictions.RMSE-mean(predictions.RMSE));
            idx = find(err==min(err),1);
    end
end
y = predictions.Y{idx};
yPred = predictions.YPred{idx};
x = 0:numel(y)-1;
plot(x,y,x,yPred)
legend("True RUL","Predicted RUL")
xlabel("Time stamp (Test data sequence)")
ylabel("RUL (Cycles)")

title("RUL for Test engine #"+idx+ " ("+lambdaCase+" case)")
end

Ссылки