基于卷积神经网络的剩余使用寿命估计

这个例子使用了:

打开实时脚本

这个例子展示了如何使用深度卷积神经网络(CNN)预测发动机的剩余使用寿命(RUL)。[1]．深度学习方法的优点是不需要手动提取特征或对模型进行特征选择来预测RUL。此外，开发基于深度学习的RUL预测模型并不需要机器健康预测和信号处理的先验知识。

下载数据集

本例使用涡扇发动机退化模拟数据集(C-MAPSS)[２]．该zip文件包含了四套不同机组(FD001、FD002、FD003、FD004)在不同运行条件和故障模式组合下的运行到故障时间序列数据。

本例仅使用FD001数据集，该数据集进一步分为训练子集和测试子集。训练子集包含100个引擎的模拟时间序列数据。每个发动机都有几个传感器，这些传感器的值在一个连续的过程中被记录在给定的实例中。因此，记录数据的顺序在长度上有所不同，对应于一个完整的运行到故障(RTF)实例。测试子集包含100个部分序列和每个序列结束时剩余使用寿命的对应值。

下载涡扇发动机退化模拟数据集到名为“CMAPSSData.zip”的文件中，并将其解压缩到名为“CMAPSSData.zip”的文件夹中“数据”在当前目录。

文件名=“CMAPSSData.zip”；如果~存在(文件名,“文件”) url =“https://ti.arc.nasa.gov/c/6/”；websave(文件名、url);结束dataFolder =“数据”；如果~存在(dataFolder“dir”mkdir (dataFolder);结束解压缩(文件名,dataFolder)

数据文件夹现在包含有26列数字的文本文件，用空格分隔。每一行是在单个操作周期中所采集的数据的快照，每列表示一个不同的变量:

第一栏:单元号
第2列:时间戳
第3-5列:操作设置
列6-26:传感器测量1-21

预处理训练数据

使用函数加载数据localLoadData．该函数从filenamePredictors并返回一个包含训练预测器和相应响应(即RUL)序列的表。每一行代表一个不同的引擎。

fileametrainpredictors = fullfile(数据文件夹，“train_FD001.txt”）;rawTrain = localLoadData(fileametrainpredictors);

检查其中一个引擎的运行到故障数据。

头(rawTrain.X {1}, 8)

ans =表8×26id时间戳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  __ _________ ____________ ____________ ____________ ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ __________________ 100 1 -0.0007 -0.0004 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 39.06 - 23.419 392 2388 100 100 2 0.0019 -0.0003 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 23.424 8131.5 8.4318 0.03 392 2388 100 39 1588 100 518.67 - 642.35 -0.0043 - 0.0003 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 518.67 642.35 1582.8 0.0007 23.344 1 4 0 1001401.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 -0.0019 - -0.0002 38.88 - 23.374 392 2388 100 1 100 554 518.67 642.37 1582.8 1406.2 14.62 21.61 2388.1 9055.1 1.3 47.28 - 522.19 2388 8133.8 8.4294 0.03 393 2388 100 38.9 23.404 -0.0043 - -0.0001 100 518.67 642.1 1584.5 1398.4 14.62 21.61 554.67 2388 2388 9049.7 1.3 47.16 521.68 8132.9 8.4108 - 0.03 38.98 - 23.367 391 2388 100 1 100 0.001 0.0001 518.67 642.48 1592.3 1397.8 9059.1 1.3 47.36 14.62 21.61 554.34 2388522.32 2388 8132.3 8.3974 0.03 392 2388 100 39.1 23.377 18 -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×1191 190 189 188 187 186 185 184

可视化一些预测器的时间序列数据。

stackedplot (rawTrain。X{1},(3、5、6、7、8、15、16、24),“XVariable”，“时间戳”）

删除变异性较小的特性

在所有时间步骤中保持不变的特征会对训练产生负面影响。使用prognosability函数来测量故障时特征的可变性。

prog =可预测性(rawTrain. prog =可预测性)X,“时间戳”）;

可以观察到，对于某些特征，可预测性等于零或NaN。放弃这些特性。

idxToRemove = prog。| isnan(program .Variables);featToRetain = prog.Properties.VariableNames(~idxToRemove);为i = 1:height(rawTrain) rawTrain。X{i} = rawTrain.X{i}{:，featToRetain};结束

规范化训练预测器

将训练预测器归一化，使其具有零均值和单位方差。

[~，Xmu,Xsigma] = zscore(vertcat(rawTrain.X{:}));preTrain = table();为i = 1:数字(rawTrain.X)预训练。X{i} = (rawTrain。X{i} - Xmu) ./ Xsigma;结束

剪辑的反应

该步骤是可选的。为了使网络专注于引擎更有可能故障的数据部分(引擎生命周期结束)，将响应剪辑到阈值150。这使得网络平等地对待具有更高RUL值的实例。

clipResponses = true;如果clipResponses rulThreshold = 150;为i = 1:数字(rawTrain.Y)预训练。Y｛i} = min(rawTrain.Y{i},rulThreshold);结束结束

该图显示了第一次观测和相应的剪切响应。

为CNN培训准备数据

本例中使用的CNN架构以2D格式(类似于图像)接收输入数据，其中一个维度表示序列长度，另一个维度表示特征的数量。然而，每个引擎的时间序列数据在长度上是不同的。因此，从时间序列数据中提取固定长度的序列，使其与CNN训练兼容。

生成固定长度的序列

使用localGenerateSequences函数从每个引擎的时间序列数据生成固定窗口长度的序列。每个序列的对应响应变量表示该序列的最后一个时间戳的RUL。

WindowLength描述在给定时间戳中用于预测RUL的时间序列数据的长度
步测量的是两个连续窗口长度之间的位移吗WindowLength

为了更好地理解localGenerateSequences工作，考虑一个长度为10的序列数据。如果所选参数WindowLength4和步的1被用来创建序列localGenerateSequences会在一个表中返回7个长度为4的序列．

WindowLength40和步其中1是根据实验选择的。它们可以设置为其他值，以研究改进的模型性能。

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

重塑imageInputLayer的序列数据

的imageInputLayer(深度学习工具箱)在CNN体系结构中使用的方法需要输入数据的大小，指定为整数行向量[h w c],在那里h，w,c分别对应通道的高度、宽度和数量。因此，将数据重塑为具有3个维度[WindowLength numFeatures 1]．

numFeatures = size(dataTrain.X{1}，2);InputSize = [WindowLength numFeatures 1];为i = 1:size(dataTrain,1) dataTrain。X{i} =重塑(dataTrain.X{i}，InputSize);结束

网络体系结构

中描述了用于RUL估计的深度卷积神经网络结构[1]．

输入数据以2D格式准备，第一个维度表示时间序列的长度，第二个维度表示所选特征的数量。卷积层堆叠用于特征提取，然后是全连接层。

所选网络结构仅沿时间序列方向应用1D卷积。这意味着特征的顺序InputSize不会影响训练，每次只考虑一个特征的趋势。

定义网络架构。创建一个CNN，由CNN层和relu层的四个连续集组成filterSize而且numFilters的两个输入参数convolution2dLayer，然后是一个完全连接的大小层numHiddenUnits和一个退出概率为0.5的退出层。的第二个维度的值filterSize为1，则得到一维卷积。由于网络预测的是涡扇发动机的剩余使用寿命(RUL)，所以设置numResponses为1。

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

培训网络

指定trainingOptions(深度学习工具箱)．使用“adam”求解器，用512大小的小批量训练20个epoch。集LearnRateSchedule来分段降低训练期间的学习率LearnRateDropFactor到0.3作为训练的下降因子。指定学习率0.003以减慢学习速度，并在每个epoch对数据进行洗牌，使网络健壮。打开训练进度图，并关闭命令窗口输出(详细的）.

maxEpochs = 20;miniBatchSize = 512;选项= trainingOptions(“亚当”，.．.“LearnRateSchedule”，“分段”，.．.“LearnRateDropFactor”, 0.3,.．.“MaxEpochs”maxEpochs,.．.“MiniBatchSize”miniBatchSize,.．.“InitialLearnRate”, 0.003,.．.“洗牌”，“every-epoch”，.．.“阴谋”，“训练进步”，.．.“详细”, 0);

使用trainNetwork．

net = trainNetwork(dataTrain,layers,options);

绘制网络的层图，以可视化底层网络架构。

图;lgraph = layerGraph(net.Layers);情节(lgraph)

测试网络

测试数据包含100个部分序列，以及每个序列末尾的剩余使用寿命对应值。

filenameTestPredictors = fullfile(数据文件夹，“test_FD001.txt”）;filenameTestResponses = fullfile(数据文件夹，“RUL_FD001.txt”）;dataTest = localLoadData(filenameTestPredictors,filenameTestResponses);

通过执行与训练数据集相同的预处理步骤，为预测准备测试数据集。

为i = 1:numel(dataTest. x) dataTest. xX{i} = dataTest.X{i}{:，featToRetain};人数(。X{i} = (dataTest。X{i} - Xmu) ./ Xsigma;如果clipResponses人数(。Y｛i} = min(dataTest.Y{i},rulThreshold);结束结束

该网络将无法对长度小于的时间序列进行预测WindowLength．因此，删除序列小于的测试单元WindowLength．

unitlength = 0(数字(dataTest.Y)，1);为i = 1:数字(dataTest.Y) unitlength (i) =数字(dataTest.Y{i，:});结束dataTest(unitlength

创建一个表来存储预测的响应(YPred)连同真实的回应(Y）.

预测=表(“大小”(高度(人数()2),“VariableTypes”，[“细胞”，“细胞”),“VariableNames”，[“Y”，“YPred”]);为i=1:height(dataTest) unit = localGenerateSequences(dataTest(i，:)，WindowLength,Stride);预测。Y｛i} = unit.Y; predictions.YPred{i} = predict(net,unit,“MiniBatchSize”, miniBatchSize);结束

性能指标

计算测试序列的所有时间周期的均方根误差(RMSE)，以比较网络在测试数据上的表现。

为i = 1:大小(预测，1)预测。rmse (i) =平方根(均值(预测。Y{i} - predictions.YPred{i}).^2));结束

下面的直方图有助于可视化RMSE值在所有测试引擎中的分布。

图;柱状图(预测。RMSE,“NumBins”10);标题(RMSE(均数:+ round(mean(prediction . rmse)，2) +， StDev: "+ round(std(prediction . rmse)，2) +“)”）;ylabel (“频率”）;包含(“RMSE”）;

此外，看看网络预测器在测试引擎中给定的数据序列中是如何执行的。使用localLambdaPlot函数绘制预测RUL相对于任何测试引擎的真实RUL。

图;localLambdaPlot(预测,“随机”）;

结果表明，CNN深度学习结构用于估计涡轮发动机数据的RUL是一种替代方法。所有时间戳上的RMSE值表明网络能够在给定测试序列数据结束时表现良好。这表明，在试图预测RUL时，有一些传感器值的小历史是很重要的。

辅助函数

数据函数

该函数从提供的文本文件加载运行到故障的数据，将时间序列数据及其对应的RUL值分组为表中的预测器和响应。

函数data = localLoadData(filenamePredictors,varargin)如果isempty(varargin) filenamerresponses = [];其他的filenameResponses = varargin{:};结束将文本文件作为表格加载rawData = readtable(fileamepredictors);向表中添加变量名VarNames = {.．.“id”，“时间戳”，“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;如果~isempty(filesameresponses) RULTest = dlmread(filesameresponses);结束为每个单元ID分割信号id = rawData{:，1};nID =唯一的(id);numObservations = numel(nID);初始化一个表来存储数据数据=表(“大小”[numObservations 2],.．.“VariableTypes”, {“细胞”，“细胞”}，.．.“VariableNames”, {“X”，“Y”})；为i=1:numObservations idx = IDs == nID(i);数据。X{i} = rawData(idx，:);如果isempty (filenameResponses)%从列车数据的时间列计算RUL数据。Y｛i} = flipud(rawData.timeStamp(idx))-1;其他的%使用来自filenamerresponses的RUL值作为测试数据sequenceLength = sum(idx);endRUL = RULTest(i);数据。Y｛i} = [endRUL+sequenceLength-1:-1:endRUL]';% #好< NBRAK >结束结束结束

生成序列函数

该函数从给定的时间序列数据生成序列WindowLength而且步．输出是一个以矩阵表示的序列表，以向量表示的对应RUL值。

函数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”、假);seqTable。X = vertcat(temp{:});temp = cellfun(getY,dataTable。Y,“UniformOutput”、假);seqTable。Y= vertcat(temp{:});结束%子函数函数seqCell = generateSequences(tsData,WindowLength,Stride)%使用WindowLength和Stride从时间序列数据中返回序列单元格数组创建一个函数来提取给定起始索引的单个序列getSeq = @(idx) tsData(1+idx:WindowLength+idx，:);确定序列的起始索引idxShift = num2cell(0:Stride:size(tsData,1)-WindowLength)';%提取序列seqCell = cellfun(getSeq,idxShift，“UniformOutput”、假);结束

Plot函数

此帮助函数接受预测表和lambdaCase，在整个序列中(在每个时间戳)绘制预测的RUL与真实的RUL，以便更好地显示预测如何随着每个时间戳而变化。第二个论点lambdaCase该函数可以是测试引擎号或一组有效字符串，以查找引擎号，如“随机”、“最佳”、“最差”或“平均”。

函数lambdaCase localLambdaPlot(预测)如果isnumeric(lambdaCase) idx = lambdaCase;其他的开关lambdaCase情况下｛“随机”，“随机”，“r”} idx = randperm(高度(预测)，1);随机选择一个测试用例来绘制情况下｛“最佳”，“最佳”，“b”} idx = find(预测。RMSE == min(prediction .RMSE));%最好的情况情况下｛“坏的”，“坏的”，“w”} idx = find(预测。RMSE == max(forecasting .RMSE));%最坏情况情况下｛“平均”，“平均”，“一个”} err = abs(forecasts . rmse -mean(forecasts . rmse));Idx = find(err==min(err)，1);结束结束y = predictions.Y{idx};yPred = predictions.YPred{idx};X = 0:数字(y)-1;情节(x, y, x, yPred)传奇(“真正的原则”，“预测原则”)包含(“时间戳(测试数据序列)”) ylabel (”原则(周期)”)标题(“测试引擎#规则”+ idx +”(“+ lambdaCase +“)”）结束

参考文献

李欣，丁琪，j - q。Sun，“使用深度卷积神经网络在预后中的剩余有用寿命估计，”可靠性工程与系统安全， 2018年4月，第172卷，第1-11页
Saxena, Abhinav, Kai Goebel。涡扇发动机退化仿真数据集NASA艾姆斯预测数据仓库https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/#turbofan美国宇航局艾姆斯研究中心，莫菲特场，加利福尼亚州