MATLAB环境下使用卷积神经网络对涡轮风扇发动机剩余使用寿命估计

之前基于python环境做了一些剩余使用寿命RUL估计的东西

NASA涡轮喷气发动机风扇的剩余寿命RUL预测-基于传统机器学习 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/528098659

NASA涡轮喷气发动机风扇的剩余寿命RUL预测-基于LSTM网络,Lookback=20 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/528324129

NASA涡轮喷气发动机风扇的剩余寿命RUL预测-几种不同方法的对比 - 哥廷根数学学派的文章 - 知乎
https://zhuanlan.zhihu.com/p/528451358

本篇使用卷积神经网络CNN进行涡轮风扇发动机的剩余使用寿命RUL估计，主要讲解如何使用卷积神经网络 CNN)来预测涡轮风扇发动机的剩余使用寿命 (RUL)。 深度学习的优势在于模型无需通过手动特征提取或特征选择来预测 RUL。 此外，基于深度学习的 RUL 预测模型不需要信号处理的先验知识。

本篇使用涡轮风扇发动机退化仿真数据集，包含在运行条件和故障模式的不同组合下模拟的四个不同组（即 FD001、FD002、FD003、FD004）的运行到故障时间序列数据。本例仅使用 FD001 数据集，该数据集进一步分为训练集和测试集。 训练集包含 100 个发动机引擎的模拟时间序列数据，每个引擎都有几个传感器，并且对应于完整的run-to-failure样本。 测试集包含部分序列和每个序列结束时的剩余使用寿命的对应值。

数据文件夹现在包含由空格分隔的 26 列数字的文本文件， 每一行都是在单个运行周期中获取的数据，每一列代表一个不同的变量。

Column 1: Unit number

Column 2: Time-stamp

Columns 3–5: Operational settings

Columns 6–26: Sensor measurements 1–21

数据预处理

加载数据，返回一个包含训练预测变量和相应响应（即 RUL）序列的表， 每行代表一个发动机引擎

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

检查其中一个引擎的run-to-failure数据

head(rawTrain.X{1},8)

查看其中一个引擎的响应数据

rawTrain.Y{1}(1:8)

可视化部分时间序列

删除Less Variability的特征

在所有时间步长上“保持近似不变”的特征会对训练产生负面影响。 使用prognosability函数来衡量失效时特征的variability（我也不想中文夹杂英文的）。

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

此外为了让卷积神经网络专注于引擎更可能发生故障（引擎生命周期结束）的数据部分，将响应进行clip剪裁

clipResponses = true;

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

为了最大限度地减少添加到小批量的填充量，按序列长度对训练数据进行排序。 然后，选择一个 mini-batch 大小，它可以均匀地划分训练数据并减少 mini-batch 中的填充量，按序列长度对训练数据进行排序。

for i = 1:size(preTrain,1)
    preTrain.X{i} = preTrain.X{i}';    %转置训练数据以在第一维中具有特征
    preTrain.Y{i} = preTrain.Y{i}';    %转置对应于训练数据的响应
    sequence = preTrain.X{i};
    sequenceLengths(i) = size(sequence,2); 
end

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

网络结构

输入数据以序列格式进行处理和排序，第1维表示特征数量，第2维表示时间序列长度。卷积层与批归一化层在一起，然后是激活层（relu），然后堆叠在一起进行特征提取。最后使用全连接层和回归层得到最终的 RUL 值作为输出。 所选网络架构仅沿时间序列方向应用一维（1D）卷积，这意味着特征的顺序不会影响训练，并且一次只考虑一个特征的趋势。

定义网络架构，创建一个 CNN，由5个连续的1D-卷积 、批量归一化和一个 relu 层组成，其中 filterSize 和 numFilters 作为卷积 1D-Layer 的前两个输入参数，然后是一个大小为 numHiddenUnits 的全连接层和一个 dropout 层(概率为 0.5)。由于网络预测涡扇发动机的剩余使用寿命（RUL），因此将第2个全连接层中的 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()];

训练网络

指定训练参数。 使用“adam”优化器，指定学习率 0.01。 为防止梯度爆炸，将梯度阈值设置为 1。

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 个部分序列和每个序列结束时剩余使用寿命RUL的对应值

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函数对测试数据进行预测

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

绘制预测的 RUL 与真实的 RUL

完整的代码及数据见如下链接

https://mianbaoduo.com/o/bread/mbd-Y5iXm5hr

参考文献

X. Li, Q. Ding, and J.-Q. Sun, “Remaining useful life estimation in prognostics using deep convolution neural networks,” Reliability Engineering & System Safety, vol. 172, pp. 1–11, Apr. 2018

Saxena, Abhinav, Kai Goebel. "Turbofan Engine Degradation Simulation Data Set." NASA Ames Prognostics Data Repository https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/#turbofan, NASA Ames Research Center, Moffett Field, CA