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

Этот пример использует Турбовентиляторный набор данных симуляции ухудшения 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

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

Подготовка данных для обучения CNN

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

Сгенерируйте последовательности фиксированной длины

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

Понять лучше как localGenerateSequences работает, рассмотрите данные о последовательности длины 10. Если выбранные параметры WindowLength из 4 и Stride из 1 используются, чтобы создать последовательности затем localGenerateSequences возвратил бы в общей сложности 7 последовательностей длины 4 в table.

WindowLength из 40 и Stride из 1 был выбран на основе экспериментирования. Они могут собираться в другие значения исследовать улучшенную производительность модели.

WindowLength = 40;
Stride = 1;
dataTrain = localGenerateSequences(preTrain,WindowLength,Stride);

Измените Данные о Последовательности для imageInputLayer

imageInputLayer (Deep Learning Toolbox), используемый в архитектуре CNN, ожидает размер входных данных в виде вектора-строки из целых чисел [h w c], где hW, и c соответствуйте высоте, ширине и количеству каналов соответственно. Поэтому измените данные, чтобы иметь 3 размерности как [WindowLength numFeatures 1].

numFeatures = size(dataTrain.X{1},2);
InputSize = [WindowLength numFeatures 1];
for i = 1:size(dataTrain,1)
    dataTrain.X{i} = reshape(dataTrain.X{i},InputSize);
end

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

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

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

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

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

filterSize = [5, 1];
numHiddenUnits = 100;
numResponses = 1;

layers = [
    imageInputLayer(InputSize)
    convolution2dLayer(filterSize,10)
    reluLayer()
    convolution2dLayer(filterSize,20)
    reluLayer()
    convolution2dLayer(filterSize,10)
    reluLayer()
    convolution2dLayer([3 1],5)
    reluLayer()
    fullyConnectedLayer(numHiddenUnits)
    reluLayer()
    dropoutLayer(0.5)
    fullyConnectedLayer(numResponses)
    regressionLayer()];

Обучите сеть

Задайте trainingOptions (Deep Learning Toolbox). Обучите в течение 20 эпох с мини-пакетами размера 512 использований ‘adam’ решателя. Установите LearnRateSchedule to piecewise пропустить скорость обучения во время обучения и LearnRateDropFactor к 0,3 как понижающийся фактор для обучения. Задайте скорость обучения 0.003, чтобы замедлить изучение и переставить данные каждая эпоха, чтобы сделать сеть устойчивой. Включите график процесса обучения и выключите командное окно выход (Verbose).

maxEpochs = 20;
miniBatchSize = 512;

options = trainingOptions('adam', ...
    'LearnRateSchedule','piecewise',...
    'LearnRateDropFactor',0.3,...
    'MaxEpochs',maxEpochs, ...
    'MiniBatchSize',miniBatchSize, ...
    'InitialLearnRate',0.003, ...
    'Shuffle','every-epoch',...
    'Plots','training-progress',...
    'Verbose',0);

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

net = trainNetwork(dataTrain,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

Сеть не сможет сделать предсказания для timeseries, длина которого меньше WindowLength. Поэтому удалите тестовые модули с последовательностями, меньшими, чем WindowLength.

unitLengths = zeros(numel(dataTest.Y),1);
for i = 1:numel(dataTest.Y)
    unitLengths(i) = numel(dataTest.Y{i,:});
end
dataTest(unitLengths<WindowLength+1,:) = [];

Составьте таблицу для хранения предсказанного ответа (YPred) наряду с истинным ответом (Y).

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

for i=1:height(dataTest)
    unit = localGenerateSequences(dataTest(i,:),WindowLength,Stride);
    predictions.Y{i} = unit.Y;
    predictions.YPred{i} = predict(net,unit,'MiniBatchSize',miniBatchSize);
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

Generate Sequences function

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

function seqTable = localGenerateSequences(dataTable,WindowLength,Stride)

getX = @(X) generateSequences(X,WindowLength,Stride);
getY = @(Y) Y(WindowLength:Stride:numel(Y));

seqTable = table;
temp = cellfun(getX,dataTable.X,'UniformOutput',false);
seqTable.X = vertcat(temp{:});
temp = cellfun(getY,dataTable.Y,'UniformOutput',false);
seqTable.Y = vertcat(temp{:});
end

% sub-function
function seqCell = generateSequences(tsData,WindowLength,Stride)
% returns a cell array of sequences from time-series data using WindowLength and Stride

% create a function to extract a single sequence given start index
getSeq = @(idx) tsData(1+idx:WindowLength+idx,:);
% determine starting indices for sequences
idxShift = num2cell(0:Stride:size(tsData,1)-WindowLength)';
% extract sequences
seqCell = cellfun(getSeq,idxShift,'UniformOutput',false);
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

Ссылки