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

Этот пример использует Набор данных симуляции деградации 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 и возвращает таблицу, которая содержит обучающие предикторы и соответствующие ответные (то есть 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");

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

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

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

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

The 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 (Deep Learning Toolbox). Обучите на 20 эпох с мини-батчами размера 512 с помощью решателя «adam». Задайте 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

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

Ссылки