基于相似度的剩余使用寿命估计
该示例展示了如何构建完整的剩余使用寿命(RUL)估计工作流，包括预处理、选择可趋势特征、通过传感器融合构建健康指标、训练相似度RUL估计器以及验证预测性能的步骤。本例使用来自PHM2008挑战数据集[1]的训练数据。
数据准备
由于数据集很小，将整个退化数据加载到内存中是可行的。下载数据集并解压缩到当前目录。使用helperLoadDataHelper函数，用于加载训练文本文件并将其转换为时间表单元格数组。训练数据包含260个运行到故障的模拟。这组测量被称为“集合”。
url =“https://ssd.mathworks.com/supportfiles/nnet/data/TurbofanEngineDegradationSimulationData.zip”；websave (“TurbofanEngineDegradationSimulationData.zip”url);解压缩(“TurbofanEngineDegradationSimulationData.zip”)退化数据= helperLoadData(“train_FD002.txt”）;degradationData (1:5)
ans =5×1单元格数组{149×26表}{269×26表}{206×26表}{235×26表}{154×26表}
每个集合成员都是一个有26列的表。列包含机器ID、时间戳、3个操作条件和21个传感器测量数据。
头(degradationData {1})
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 __ ____ ____________ ____________ ____________ ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ 134.998 - 0.84 100 449.44 555.32 1358.6 1137.2 5.48 194.64 2222.7 8341.9 1.02 42.02 183.06 2387.7 8048.6 9.3461 0.02 334 2223 100 14.73 8.8071 41.998 - 0.8408 100 445 549.9 1353.2 1125.8 3.91 5.71 138.51 2211.6 8304 1.02 42.2 130.42 2387.7 8072.3 9.3774 0.02 330 2212 100 10.41 6.2665 - 1 3 60 24.999 - 0.6218 462.54 537.31 1256.8 1047.5 7.05 9.02 175.71 1915.1 8001.4 0.94 36.69 164.22 2028 7864.9 10.894 0.02 309 1915 84.93 14.08 8.6723 42.008 - 0.8416 100 445 549.51 1354 1126.4 3.91 5.71 138.462211.6 8304 1.02 41.96 130.72 2387.6 8068.7 9.3528 0.02 329 2212 100 10.59 6.4701 0.6203 1 5 25 60 462.54 537.07 1257.7 1047.9 7.05 9.03 175.05 1915.1 7993.2 0.94 36.89 164.31 2028 7861.2 10.896 0.02 309 1915 84.93 14.13 8.5286 1 6 60 25.005 - 0.6205 462.54 537.02 1266.4 1048.7 7.05 9.03 175.17 1915.2 7996.1 0.94 36.78 164.27 2028 7868.9 10.891 0.02 306 1915 84.93 14.28 8.559 42.004 - 0.8409 100 445 549.74 1347.5 1127.2 3.91 5.71 138.71 2211.6 8307.8 1.02 42.19 130.49 2387.7 8075.5 9.3753 0.02330 2212 100 10.62 6.4227 1 8 20.002 0.7002 100 491.19 607.44 1481.7 1252.4 9.35 13.65 334.41 2323.9 8709.1 1.08 44.27 315.11 2388 8049.3 9.2369 0.02 365 2324 100 24.33 14.799
将退化数据分解为训练数据集和验证数据集，用于以后的性能评估。
rng (“默认”）以确保结果是可重复的numEnsemble = length(退化数据);numFold = 5;cv = cvpartition(numEnsemble，“KFold”, numFold);trainData =降解数据(training(cv, 1));validationData =退化数据(测试(cv, 1));
指定感兴趣的变量组。
varNames = string(降解数据{1}.Properties.VariableNames);timeVariable = varNames{2};条件变量= varNames(3:5);dataVariables = varNames(6:26);
可视化集成数据的示例。
Nsample = 10;figure helperPlotEnsemble(训练数据，时间变量，…[conditionVariables(1:2) dataVariables(1:2)]， nsample)
工作状态聚类
如前一节所示，在每次运行到故障的测量中，没有显示退化过程的明显趋势。在本节和下一节中，将使用操作条件从传感器信号中提取更清晰的退化趋势。
注意，每个集成成员包含3个操作条件:“op_setting_1”、“op_setting_2”和“op_setting_3”。首先，让我们从每个单元格中提取表，并将它们连接到单个表中。
trainDataUnwrap = vertcat(trainData{:});opConditionUnwrap = trainDataUnwrap(:， cellstr(conditionVariables));
在3D散点图上可视化所有操作点。它清楚地显示了6个区域，每个区域的点都非常接近。
图helperPlotClusters (opConditionUnwrap)
让我们使用聚类技术自动定位这6个聚类。这里使用K-means算法。K-means是最流行的聚类算法之一，但它会导致局部最优。用不同的初始条件重复K-means聚类算法多次，选取代价最低的结果是一个很好的做法。在本例中，该算法运行5次，结果相同。
Opts = statset(“显示”，“最后一次”）;[clusterIndex, centers] = kmeans(table2array(opConditionUnwrap)， 6，…“距离”，“sqeuclidean”，“复制”5,“选项”、选择);
重复1,1次迭代，总距离和= 0.331378。重复2,1次迭代，总距离和= 0.331378。重复3,1次迭代，总距离和= 0.331378。重复4,1次迭代，总距离和= 0.331378。重复5,1次迭代，总距离和= 0.331378。最佳距离总和= 0.331378
可视化聚类结果和识别的聚类质心。
图helperPlotClusters(opConditionUnwrap, clusterIndex, centers)
如图所示，聚类算法成功地找到了6个工作机制。
工作机制正常化
让我们对按不同工作机制分组的度量进行标准化。首先，计算每个传感器测量的平均值和标准偏差，按上一节中确定的工作状态分组。
Centerstats = struct(“的意思是”、表()“SD”表());为v = dataVariables centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v))， clusterIndex);centerstats.SD.(char(v)) = splitapply(@std, trainDataUnwrap.(char(v))， clusterIndex);结束centerstats。的意思是
ans =6×21表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 ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ 489.05 - 604.91 2319 1502 1311.2 10.52 15.493 394.33 8785.8 1.26 45.487 371.45 2388.2 8135.5 8.66450.03 369.69 2319 100 28.527 17.115 462.54 536.86 1262.7 1050.2 7.05 9.043 1915.4 8015.6 0.93992 36.798 164.57 2028.3 78792 10.913 0.02 307.34 1915 84.93 14.263 8.5587 445 549.7 1354.4 1127.7 3.91 5.7158 138.63 2212 8328.1 1.0202 42.146 130.55 2388.1 8089.7 9.3745 0.02 331.06 2212 100 10.588 6.3516 491.19 607.56 1485.5 1253 9.35 13.657 334.5 2324 8730 1.0777 44.446 314.88 2388.2 8066.1 9.2325 0.022112 365.37 2388.1 24.455 14.674 518.67 642.67 1590.39063.3 1.3 47.538 521.41 2388.1 8142 8.4419 0.03 393.18 2388 100 38.82 23.292 449.44 555.8 1366.7 1131.5 5.48 8.0003 194.45 2223 8356.4 1.0204 41.981 183.03 2388.1 8072.5 9.3319 0.02 334.22 2223 100 14.831 8.8983
centerstats。SD
ans =6×21表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 __________ ________ ________ ________ __________ _________ ________ ________ ________ __________ _________ _________ _________ _________ _________ __________ _________ _________ __________ _________ _________ e-13 e-11 0.47582 5.8057 8.5359 3.2687 1.0346 0.0047002 0.67361 0.095286 18.819 - 0.000175520.2538 0.53621 0.098213 16.507 0.038392 6.6965e-16 1.49 00 0.14441 0.08598 4.6615e- 13 0.0042209 0.44862 0.27006 14.479 0.00087595 0.20843 0.34493 0.28592 13.427 0.043297 3.4003e-16 1.2991 0 3.5246e-12 0.11218 0.066392 0 0.6464 1.1014e-13 0.0049312 0.43972 0.30478 18.325 0.0015563 0.23965 0.34011 0.3286 1.773 0.037476 1.0652e-15 1.4114 000 0.10832 0.07768 2.5049e-13 0.0047449 0.6046 0.12828 17.58 0.004193 0.238260.49843 0.13156 15.457 0.038886 0.004082 1.4729 00 0.13531 0.079592 1.2279e-11 0.4967 0.86985 0.07082 19.889 2.7314e-14 0.26508 0.73713 0.070545 17.105 0.037289 6.6619e-16 1.5332 00 0.17871 0.1065 1.4098e-11 0.4401 5.6013 7.4461 2.2828e-13 0.0017508 0.47564 0.28634 17.446 0.0019595 0.23078 0.37631 0.30766 15.825 0.03856e -16 1.4133 00 0.11347 0.067192
每个体制下的统计数据可以用来对训练数据进行归一化。对于每个集成成员，提取每行的工作点，计算其到每个簇中心的距离，并找到最近的簇中心。然后，对于每个传感器测量，减去平均值并除以该集群的标准偏差。如果标准偏差接近0，则将归一化传感器测量值设置为0，因为几乎恒定的传感器测量值对于剩余使用寿命估计没有用处。请参阅最后一节“Helper函数”，了解更多关于regenormalization函数的详细信息。
traindatanormalization = cellfun(@(data) regenormalization (data, centers, centerstats)，…trainData,“UniformOutput”、假);
根据工作状态将数据规范化。一些传感器测量的退化趋势现在在归一化后显示出来。
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(1:4)， nsample)
Trendability分析
现在选择最具趋势的传感器测量来构建用于预测的健康指示器。对于每个传感器测量，估计线性退化模型，并对信号的斜率进行排序。
numSensors = length(dataVariables);signalSlope = 0 (numSensors, 1);警告=警告(“关闭”）;为ct = 1:numSensors tmp = cellfun(@(tbl) tbl(:， cellstr(dataVariables(ct)))， trainDataNormalized，“UniformOutput”、假);mdl =线性退化模型();%创建模型fit (mdl、tmp);%列车模式signalSlope(ct) = mdl.Theta;结束警告(警告);
对信号斜率进行排序，选取斜率最大的8个传感器。
[~， idx] = sort(abs(signalSlope)，“下”）;sensorTrended = sort(idx(1:8))
sensorTrended =8×12 3 4 7 11 12 15 17
可视化选择的可趋势传感器测量值。
figure helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(sensorTrended(3:6))， nsample)
请注意，一些最具趋势的信号显示出积极的趋势，而另一些则显示出消极的趋势。
构建运行状况指示器
本节重点介绍如何将传感器测量数据融合为单个健康指标，并利用该指标训练基于相似性的模型。
假定所有运行到故障数据都以正常状态开始。开始时的运行状况状况被赋值为1，故障时的运行状况被赋值为0。假定运行状况随时间从1线性下降到0。这种线性退化用于帮助融合传感器值。文献[2-5]描述了更复杂的传感器融合技术。
为j=1:numel(trainDataNormalized) data = trainDataNormalized{j};Rul = max(data.time)-data.time;数据。Health_condition = rul / max(rul);trainDataNormalized{j} = data;结束
想象健康状况。
figure helperPlotEnsemble(trainDataNormalized, timeVariable，“health_condition”nsample)
随着降级速度的不同，所有集成成员的健康状况从1变为0。
现在拟合一个健康状况的线性回归模型，以最具趋势的传感器测量作为回归量:
健康状况~ 1 +传感器2 +传感器3 +传感器4 +传感器7 +传感器11 +传感器12 +传感器15 +传感器17
trainDataNormalizedUnwrap = vertcat(trainDataNormalized{:});sensorToFuse = dataVariables(sensorTrended);X = trainDataNormalizedUnwrap{:， cellstr(sensorToFuse)};y = trainDataNormalizedUnwrap.health_condition;regModel = fitlm(X,y);偏差= regModel.Coefficients.Estimate(1)
偏差= 0.5000
权重= regModel.Coefficients.Estimate(2:end)
重量=8×1-0.0323 -0.0300 -0.0527 0.0057 -0.0646 0.0054 -0.0431 -0.0379
通过将传感器测量值与其相关权重相乘来构造单个运行状况指示器。
trainDataFused = cellfun(@(data)降解sensorfusion (data, sensorToFuse, weights)， trainDataNormalized，…“UniformOutput”、假);
可视化训练数据的融合运行状况指示器。
图helperPlotEnsemble(trainDataFused， []， 1, nsample) xlabel(“时间”) ylabel (“健康指示器”)标题(的训练数据）
来自多个传感器的数据被融合成一个单一的健康指示器。运行状况指示器通过移动平均过滤器平滑。更多细节请参见最后一节中的辅助函数“降解sensorfusion”。
对验证数据应用相同的操作
对验证数据集重复状态归一化和传感器融合过程。
validationDataNormalized = cellfun(@(data) regenormalization (data, centers, centerstats)，…validationData,“UniformOutput”、假);validationDataFused = cellfun(@(数据)降解sensorfusion(数据，sensorToFuse，权重)，…validationDataNormalized,“UniformOutput”、假);
查看验证数据的运行状况指示器。
图helperPlotEnsemble(validationDataFused， []， 1, nsample) xlabel(“时间”) ylabel (“健康指示器”)标题(验证数据的）
建立相似RUL模型
现在利用训练数据建立一个基于残差的相似度RUL模型。在这种情况下，模型尝试用二阶多项式拟合每个融合数据。
数据之间的距离 和数据 是由残差的1范数计算的吗
在哪里 是否是机器的运行状况指示器 ， 机器的估计运行状况指标是 利用二阶多项式模型在机器中辨识 ．
相似度分数由下式计算
在验证数据集中给定一个集合成员，模型将在训练数据集中找到最近的50个集合成员，根据这50个集合成员拟合一个概率分布，并使用分布的中位数作为RUL的估计。
mdl = residualSimilarityModel(…“方法”，“poly2”，…“距离”，“绝对”，…“NumNearestNeighbors”, 50岁,…“标准化”1);fit (mdl trainDataFused);
绩效评估
为了评估相似RUL模型，使用样本验证数据的50%、70%和90%来预测其RUL。
断点= [0.5,0.7,0.9];validationdatafmp = validationDataFused{3};%使用一个验证数据进行说明
使用第一个断点之前的验证数据，即生命周期的50%。
Bpidx = 1;validationDataTmp50 = validationDataTmp(1:ceil(end*断点(bpidx))，:);trueRUL = length(validationDataTmp) - length(validationDataTmp50);[estRUL, ciRUL, pdfRUL] = predictRUL(mdl, validationDataTmp50);
可视化截断50%的验证数据及其最近的邻居。
figure compare(mdl, validationDataTmp50);
将估计的RUL与真实的RUL以及估计的RUL的概率分布进行比较。
图helpperplotruldistribution (trueRUL, estRUL, pdfRUL, ciRUL)
当机器处于中等健康状态时，预估RUL与真实RUL之间存在较大误差。在本例中，最相似的10条曲线在开始时接近，但在接近失效状态时分叉，导致RUL分布中大致有两种模式。
使用第二个断点之前的验证数据，即生命周期的70%。
Bpidx = 2;validationDataTmp70 = validationDataTmp(1:ceil(end*断点(bpidx))，:);trueRUL = length(validationDataTmp) - length(validationDataTmp70);[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp70);figure compare(mdl, validationDataTmp70);
图helpperplotruldistribution (trueRUL, estRUL, pdfRUL, ciRUL)
当观察到更多的数据时，RUL估计会得到增强。
使用第三个断点之前的验证数据，即生命周期的90%。
Bpidx = 3;validationDataTmp90 = validationDataTmp(1:ceil(end*断点(bpidx))，:);trueRUL = length(validationDataTmp) - length(validationDataTmp90);[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp90);figure compare(mdl, validationDataTmp90);
图helpperplotruldistribution (trueRUL, estRUL, pdfRUL, ciRUL)
当机器接近故障时，RUL估计在本例中得到了更强的增强。
现在，对整个验证数据集重复相同的评估过程，并计算每个断点的估计RUL与真实RUL之间的误差。
numValidation = length(validationDataFused);numBreakpoint =长度(断点);错误=零(numValidation, numBreakpoint);为dataIdx = 1:numValidation tmpData = validationDataFused{dataIdx};为bpidx = 1:numBreakpoint tmpDataTest = tmpData(1:ceil(end*breakpoint(bpidx))，:);trueRUL = length(tmpData) - length(tmpDataTest);[estRUL， ~， ~] = predictRUL(mdl, tmpDataTest);error(dataIdx, bpidx) = estRUL - trueRUL;结束结束
可视化每个断点的误差直方图及其概率分布。
[pdf50, x50] = ksdensity(error(:， 1));[pdf70, x70] = ksdensity(error(:， 2));[pdf90, x90] = ksdensity(error(:， 3));图ax(1) = subplot(3,1,1);持有在直方图(错误(:1),“BinWidth”5,“归一化”，“pdf”) plot(x50, pdf50)等待从包含(预测误差的)标题(“使用每个验证集成成员的前50%的RUL预测错误”) ax(2) = subplot(3,1,2);持有在直方图(错误(:,2),“BinWidth”5,“归一化”，“pdf”) plot(x70, pdf70) hold住从包含(预测误差的)标题(“使用每个验证集成成员的前70%的RUL预测错误”) ax(3) = subplot(3,1,3);持有在直方图(错误(:,3),“BinWidth”5,“归一化”，“pdf”) plot(x90, pdf90)等待从包含(预测误差的)标题(“使用每个验证集成成员的前90%的RUL预测错误”) linkaxes (ax)
在箱形图中绘制预测误差，以可视化中位数、25-75分位数和异常值。
图箱线图(错误,“标签”, {“50%”，“70%”，“90%”}) ylabel (预测误差的)标题(“使用每个验证集合成员的不同百分比的预测错误”）
计算和可视化预测误差的平均值和标准偏差。
errorMean = mean(错误)
errorMean =1×3-1.1217 9.5186 9.6540
errorMedian =中位数(错误)
errorMedian =1×31.3798 11.8172 10.3457
errorSD = std(错误)
errorSD =1×321.7315 13.5194 12.3083
figure errorbar([50 70 90]， errorMean, errorSD，“o”，“MarkerEdgeColor”，“r”xlim([40,100]) xlabel(“用于RUL预测的验证数据的百分比”) ylabel (预测误差的)传说(“1个标准差误差条的平均预测误差”，“位置”，“南”）
结果表明，当观测到的数据越多，误差就越集中在0附近(离群值越少)。
参考文献
[1] A. Saxena和K. Goebel(2008)。“PHM08挑战数据集”，NASA艾姆斯预测数据库(http://ti.arc.nasa.gov/project/prognostic-data-repository)， NASA艾姆斯研究中心，莫菲特场，加州
[2]罗默尔，迈克尔J.，格雷戈里J.卡普钦斯基，迈克尔H.舍勒。“利用健康管理信息融合改进诊断和预后评估。”AUTOTESTCON程序，2001。IEEE系统准备技术会议．IEEE 2001。
格贝尔，凯，皮耶罗·博尼松。“恒负荷系统的预测信息融合。”信息融合，2005年第8届国际学术会议．卷。2。IEEE 2005。
[4]王，彭，大卫W.科伊特。“基于退化建模的具有多个退化措施的系统可靠性预测。”可靠性和可维护性，2004年度专题讨论会．IEEE 2004。
[5] Jardine, Andrew KS，林大明和Dragan Banjevic。机械诊断和预测实施基于状态的维护综述机械系统和信号处理20.7(2006): 1483-1510。
辅助函数
函数data = regenormalization (data, centers, centerstats)对数据的每个观察(行)进行归一化%根据观测值所属的聚类。conditionIdx = 3:5;dataIdx = 6:26;%按行操作data{:， dataIdx} = table2array(…rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats)，…数据,“SeparateInputs”、假));结束函数rowNormalized = localNormalize(row, conditionIdx, dataIdx, centers, centerstats)%每一行的归一化操作。获取工作点和传感器测量值ops = row(1, conditionIdx);sensor = row(1, dataIdx);找出离我们最近的群集中心。Dist = sum((centers - ops)。2) ^ 2;[~， idx] = min(dist);用聚类的均值和标准差归一化传感器测量值。将NaN和Inf重新分配为0。rowNormalized = (sensor - centerstats.)Mean{idx，:}) ./ centerstats。SD {idx:};rowNormalized(isnan(rowNormalized) | isinf(rowNormalized)) = 0;结束函数dataFused =降解sensorfusion(数据，sensorToFuse，权重)结合来自不同传感器的测量结果%为权重，平滑融合数据，偏移数据%，以便所有数据都从1开始根据权重融合数据dataToFuse = data{:， cellstr(sensorToFuse)};dataFusedRaw = datafususe *weights;%用移动均值平滑融合数据stepBackward = 10;stepForward = 10;dataFused = movmean(dataFusedRaw， [stepBackward stepForward]);将数据偏移为1dataFused = dataFused + 1 - dataFused(1);结束
另请参阅
相关的话题
您也可以从以下列表中选择一个网站: