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

В этом примере используется Набор данных моделирования деградации турбовентиляторного двигателя (C-MAPSS) [2]. ZIP-файл содержит данные временных рядов от запуска до отказа для четырех различных наборов (а именно FD001, FD002, FD003, FD004), моделируемых при различных комбинациях рабочих условий и режимов отказов.

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

Загрузите набор данных моделирования деградации турбовентилятора в файл с именем «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 и возвращает таблицу, которая содержит обучающие предикторы и соответствующие последовательности ответа (т.е. 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

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

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

Удаление элементов с меньшей вариабельностью

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

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

Можно наблюдать, что для некоторых признаков прогностическая способность равна нулю или 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 формате (аналогично изображению), где одно измерение представляет длину последовательности, а другое измерение представляет количество признаков. Однако временные ряды данных для каждого механизма различаются по длине. Поэтому последовательности фиксированной длины извлекаются из данных временных рядов, чтобы сделать их совместимыми для обучения CNN.

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

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

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

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

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

Изменение формы данных последовательности для imageInputLayer

imageInputLayer (Deep Learning Toolbox), используемый в архитектуре CNN, ожидает размера входных данных, указанных как вектор строк целых чисел [h w c], где h, w, и 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 (инструментарий глубокого обучения). Тренируйтесь на 20 эпох с мини-атчами размером 512 с помощью решателя «адам». Набор LearnRateSchedule кому 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

Сеть не сможет делать прогнозы для временных рядов, длина которых меньше 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

Эта функция загружает данные от запуска до отказа из предоставленного текстового файла, группирует данные временных рядов и соответствующие им значения 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

Эта функция генерирует последовательности на основе данных временных рядов с учетом 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

Ссылки