畳み込みニューラルネットワークを使用した残存耐用期間の推定

この例では次を使用します。

この例では、深層畳み込みニューラルネットワーク (CNN) [1] を使用してエンジンの残存耐用期間 (RUL) を予測する方法を説明します。深層学習の手法の利点は、モデルが RUL を予測するにあたって、特徴抽出や特徴選択を手動で行う必要がないことです。さらに、機械の健全性の予知または信号処理の事前知識は、深層学習に基づく RUL 予測モデルの構築に必須ではありません。

データセットのダウンロード

この例では、Turbofan Engine Degradation Simulation データセット [1] を使用します。データセットは zip ファイル形式であり、4 つの異なるセットについて故障に至るまで実行された時系列データ (FD001、FD002、FD003、FD004) が含まれており、動作条件と故障モードの異なる組み合わせでシミュレートされています。

この例では、FD001 データセットのみを使用しますが、これはさらに学習用とテスト用のサブセットに分割されます。学習用サブセットにはエンジン 100 台のシミュレートされた時系列データが含まれています。各エンジンにはいくつかのセンサーが搭載され、それらの値は連続プロセス内の指定されたタイミングで記録されています。そのため、記録された一連のデータの長さはさまざまで、故障に至るまで実行された (RTF) インスタンス全体に対応します。テスト用サブセットには 100 個の部分シーケンスと、各シーケンスの最後における対応する残存耐用期間の値が含まれています。

Turbofan Engine Degradation Simulation データセットを保存するディレクトリを作成します。

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

Turbofan Engine Degradation Simulation データセットをダウンロードして抽出します。

filename = matlab.internal.examples.downloadSupportFile("nnet","data/TurbofanEngineDegradationSimulationData.zip");
unzip(filename,dataFolder)

データフォルダーには、これで、スペースで区切られた 26 個の数値列からなるテキストファイルが含まれています。各行は単一の動作サイクルで取得されたデータのスナップショットであり、各列は別々の変数を表します。

列 1 — ユニット番号
列 2 — タイムスタンプ
列 3 ～ 5 — 動作設定
列 6 ～ 26 — センサー測定値 1 ～ 21

学習データの前処理

関数 localLoadData を使用してデータを読み込みます。この関数は、データファイルからデータを抽出して、学習予測子と、それに対応する応答 (RUL) のシーケンスが格納された table を返します。各行は異なるエンジンを表します。

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

1 つのエンジンの故障に至るまで実行されたデータを確認します。

head(rawTrain.X{1},8)

    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

1 つのエンジンの応答データを確認します。

rawTrain.Y{1}(1:8)

いくつかの予測子の時系列データを可視化します。

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

Figure contains an object of type stackedplot.

変動性が小さい特徴の削除

すべてのタイムステップで一定に留まっている特徴は、学習に悪影響を与える可能性があります。関数prognosability (Predictive Maintenance Toolbox)を使用して、故障時の特徴の変動性を測定します。

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

いくつかの特徴では、予知可能性が 0 と等しいか 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

応答のクリップ

応答データは、各エンジンの寿命に対する RUL 値を表し、個々のエンジン寿命に基づいています。シーケンスは、初期の測定時からエンジンの故障時までの線形劣化を想定しています。

エンジンが故障する可能性が高いデータ部分 (エンジンの耐用寿命末期) にネットワークを集中させるために、しきい値 150 で応答をクリップします。応答のクリップにより、ネットワークは RUL 値がこれを超えるインスタンスを等価として扱うようになります。

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

この図は、最初の観測値と、しきい値でクリップされた対応する応答 (RUL) を示しています。緑色のオーバーレイは、センサーと RUL の両方のプロットでクリップ領域を定義しています。

パディングに対するデータの準備

このネットワークは、各種シーケンス長の入力データをサポートしています。ネットワークを介してデータを渡す際、各ミニバッチ内のすべてのシーケンスが指定された長さになるように、ソフトウェアにより、シーケンスのパディング、切り捨て、分割が行われます。

ミニバッチに加えられるパディングの量を最小化するために、学習データをシーケンス長によって並べ替えます。次に、学習データを均等に分割してミニバッチでのパディングの量を減らすように、ミニバッチのサイズを選択します。

学習データをシーケンス長によって並べ替えます。

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] で説明されています。

ここで、最初の次元が選択した特徴の数、2 番目の次元が時間シーケンスの長さを表すシーケンス形式でデータを処理し、並べ替えを行います。畳み込み層をバッチ正規化層とまとめ、活性化層 (ここでは relu) を続け、特徴抽出のために層を一緒にスタックします。全結合層は、最終の RUL 値を出力として得るために、最後に使用されます。

選択したネットワークアーキテクチャでは、1 次元の畳み込みは時間シーケンスの順序に沿ってのみ適用されます。したがって、特徴の順序は学習に影響を与えず、一度に 1 つの特徴のトレンドのみが考慮されます。

ネットワークアーキテクチャを定義します。増加する filterSize と numFilters を convolution1dLayer の最初の 2 つの入力引数として、畳み込み 1-d、バッチ正規化、および relu 層の 5 つの連続したセットで構成される CNN を作成し、続いてサイズが numHiddenUnits の全結合層と、ドロップアウトの確率が 0.5 のドロップアウト層を作成します。このネットワークではターボファンエンジンの残存耐用期間 (RUL) を予測するため、2 番目の全結合層で numResponses を 1 に設定します。

学習データ内の変化する時間シーケンスを補正するために、convolution1dLayer で名前と値のペアの入力引数として Padding="causal" を使用します。

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)];

ネットワークの学習

trainingOptionsを指定します。Adam オプティマイザーを使用し、ミニバッチサイズを 16 として、学習を 40 エポック行います。LearnRateSchedule を piecewise に設定します。学習率を 0.01 として指定します。勾配爆発を回避するために、勾配のしきい値を 1 に設定します。長さによるシーケンスの並べ替えを維持するために、'Shuffle' を 'never' に設定します。学習の進行状況プロットを有効にし、Metrics を 'rmse' に設定して RMSE プロットを表示します。学習データには、最初の次元に特徴、2 番目の次元に時系列シーケンス、3 番目の次元にデータのバッチが含まれているため、InputDataFormats を CTB に設定します。コマンドウィンドウ出力 (Verbose) をオフにします。

maxEpochs = 40;
miniBatchSize = 16;

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

trainnet を使用してネットワークを学習させ、平均二乗誤差 (mse) として損失関数を指定します。約 1 ～ 2 分で終了します。

net = trainnet(XTrain,YTrain,layers,"mse",options);

ネットワークの層グラフをプロットして、背後にあるネットワークアーキテクチャを可視化します。

figure;
plot(net)

Figure contains an axes object. The axes object contains an object of type graphplot.

ネットワークのテスト

テストデータには 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;
    dataTest.Y{i} = min(dataTest.Y{i},rulThreshold);
end

予測された応答 (YPred) と真の応答 (Y) を保存するための table を作成します。minibatchpredict を使用して、テストデータについての予測を行います。関数によってパディングがテストデータに追加されるのを防ぐには、テストデータ中の各観測値の RUL 値を得られるよう、ミニバッチのサイズに 1 を指定します。テストデータも学習データと同じ形式であるため、InputDataFormats を CTB に設定します。

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} = minibatchpredict(net,unit,MiniBatchSize=1,InputDataFormats='CTB');
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');

Figure contains an axes object. The axes object with title RMSE ( Mean: 17.53 , StDev: 7.02 ), xlabel RMSE, ylabel Frequency contains an object of type histogram.

また、ネットワークの予測子がテストエンジンの特定のシーケンスデータ全体でどのような精度を示しているかを確認するには、関数 localLambdaPlot を使用して、ランダムなテストエンジンの真の RUL に対し、予測された RUL をプロットします。

figure;
localLambdaPlot(predictions,"random");

Figure contains an axes object. The axes object with title RUL for Test engine #59 (random case), xlabel Time stamp (Test data sequence), ylabel RUL (Cycles) contains 2 objects of type line. These objects represent True RUL, Predicted RUL.

この結果から、ターボエンジンデータの RUL を推定するための CNN 深層学習アーキテクチャは RUL を予測するための実行可能なアプローチであることがわかります。すべてのタイムスタンプの RMSE 値は、特定のテストシーケンスデータの終盤までネットワークの精度が高かったことを示しています。したがって、RUL の予測を試みるには、センサーデータの簡単な履歴を用意するのが重要です。

補助関数

データ関数の読み込み

この関数は、指定されたテキストファイルから故障に至るまで実行されたデータを読み込み、時系列データとそれに対応する RUL 値を、予測子と応答として table 内でグループ化します。

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 = readmatrix(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

ラムダプロット関数

この補助関数は table predictions と lambdaCase 引数を受け取り、予測された RUL を、真の RUL に対しシーケンス全体を通して (タイムスタンプごとに) プロットします。これにより、タイムスタンプごとの予測の推移が可視化されます。2 番目の引数 lambdaCase は、テストエンジンの番号、あるいは "random"、"best"、"worst"、または "average" など、エンジンの番号を見つけるための有効な string のセットの 1 つとすることができます。

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

参考文献

Li, Xiang, Qian Ding, and Jian-Qiao Sun."Remaining Useful Life Estimation in Prognostics Using Deep Convolution Neural Networks."Reliability Engineering & System Safety 172 (April 2018):1–11. https://doi.org/10.1016/j.ress.2017.11.021.

参考

imageInputLayer | prognosability (Predictive Maintenance Toolbox) | trainingOptions

畳み込みニューラルネットワークを使用した残存耐用期間の推定

データセットのダウンロード

学習データの前処理

ネットワークアーキテクチャ

ネットワークの学習

ネットワークのテスト

補助関数

データ関数の読み込み

ラムダプロット関数

参考文献

参考

関連するトピック

外部の Web サイト

畳み込みニューラル ネットワークを使用した残存耐用期間の推定

データ セットのダウンロード

学習データの前処理

ネットワーク アーキテクチャ

ネットワークの学習

ネットワークのテスト

補助関数

データ関数の読み込み

ラムダ プロット関数

参考文献

参考

関連するトピック

外部の Web サイト

畳み込みニューラルネットワークを使用した残存耐用期間の推定

データセットのダウンロード

ネットワークアーキテクチャ

ラムダプロット関数