基于相似性的剩余可用寿命估计

打开生活的脚本

该示例显示了如何构建完整的剩余使用寿命（RUL）估计工作流，包括预处理、选择趋势特征、通过传感器融合构建健康指标、培训相似性RUL估计器以及验证预测性能的步骤。该示例使用来自的PHM2008挑战数据集的训练数据https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/[1]。

数据准备

由于数据集是小它是将整个降解的数据加载到存储器可行的。下载并解压缩从数据集https://ti.arc.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/转到当前目录。使用helperLoadData辅助函数用于加载训练文本文件并将其转换为一个单元格数组的时间表。训练数据包含218个从运行到失败的模拟。这组测量值称为“集合”。

degradationData = helperLoadData (“train.txt”）;degradationData (1:5)

ans =5×1单元阵列{223×26表} {164×26表} {150×26表} {159×26表} {357×26表}

每个集合成员是一个包含26列的表。列包含机器ID，时间戳，3个操作条件和21个传感器测量的数据。

头（degradationData {1}）

ans =表8×26本次研究的目的是：传感器传感器1传感器2传感器3传感器传感器2传感器3传感器传感器3传感器传感器4传感器5传感器5传感器6传感器传感器6传感器6传感器7传感器8传感器8传感器8传感器8传感器8传感器8传感器8传感器8传感器8传感器9传感器10传感器8传感器8传感器8传感器8传感器9传感器10传感器10传感器11传感器12传感器12传感器12传感器13传感器14传感器14传感器15传感器15传感器16传感器17传感器18传感器18传感器18传感器19传感器传感器传感器19传感器20传感器传感器传感器传感器20传感器传感器传感器传感器21个传感器传感器传感器21个传感器传感器传感器传感器21个传感器传感器。这是从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从从(小标题)(小标题)(小标题)(小标题)(小标题)(小标题)(小标题)(小标题)(小标题)________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ _________ 1 1 10.005 0.2501 20 489.05 604.13 1499.5 1310 10.52 15.49 394.88 2318.9 8770.2 1.26 45.4 372.15 2388.1 8120.8 8.6216 0.03 368 2319 100 28.58 17.174 1 2 0.0015 0.0003 100 518.67 642.13 1584.5 1404 14.62 21.61 553.67 2388 9045.8 1.3 47.29 521.81 2388.2 8132.9 8.3907 0.03 391 2388 100 38.99 23.362 1 3 34.999 0.8401 60 449.44 555.42 1368.2 1122.5 5.48 8 194.93 2222.9 8343.9 1.02 41.92 183.26 2387.9 8063.8 9.3557 0.02 334 2223 100 14.83 8.8555 1 4 20.003 0.7005 0 491.19 607.03 1488.4 1249.2 9.35 13.65 334.82 2323.8 8721.5 1.08 44.26 314.84 2388.1 8052.3 9.2231 0.02 364 2324 100 24.42 14.783 1 5 42.004 0.8405 40 445 549.52 1354.5 1124.3 3.91 5.71 138.24 2211.8 8314.6 1.02 41.79 130.44 2387.9 8083.7 9.2986 0.02 330 2212 100 10.99 6.4025 1 6 20.003 0.7017 0 491.19 607.37 1480.5 1258.9 9.35 13.65 334.51 2323.9 8711.4 1.08 44.4 315.36 2388.1 8053.2 9.2276 0.02 364 2324 100 24.44 14.702 1 7 42 0.84 40 445 549.57 1354.4 1131.4 3.91 5.71 139.11 2211.8 8316.9 1.02 42.09 130.16 2387.9 8082 9.3753 0.02 331 2212 100 10.53 6.4254 1 8 0.0011 0 100 518.67 642.08 1589.5 1407.6 14.62 21.61 553.48 2388.1 9050.4 1.3 47.5 521.74 2388 8133.3 8.4339 0.03 391 2388 100 38.98 23.234

将降级数据分解为训练数据集和验证数据集，以便稍后进行性能评估。

RNG（'默认'）%以确保结果的可重复性numEnsemble =长度(degradationData);numFold = 5;简历= cvpartition (numEnsemble,“KFold”，numFold）;trainData = degradationData（训练（CV，1））;validationData = degradationData（试验（CV，1））;

指定感兴趣的变量组。

字符串(degradationData varNames = {1} .Properties.VariableNames);timeVariable = varNames {2};conditionVariables = varNames (3:5);dataVariables = varNames(第一);

可视化集成数据的示例。

nsample = 10;图helperPlotEnsemble (trainData timeVariable,．..[conditionVariables (1:2) dataVariables (1:2)], nsample)

工作制度集群

正如上一节中所示，是表示各运行到失败测量降解过程没有明显的趋势。在这个和下一个部分，操作条件将被用于提取从传感器信号更清晰的退化的趋势。

注意，每个ensemble成员包含3个运行条件:“op_setting_1”、“op_setting_2”和“op_setting_3”。首先，让我们从每个单元格中提取表，并将它们连接到单个表中。

trainDataUnwrap=vertcat（列车数据{：}）；opConditionUnwrap=trainDataUnwrap（：，cellstr（conditionVariables））；

在3D散点图上可视化所有操作点。它清楚地显示了6个区域，每个区域的点都非常接近。

图helperPlotClusters（opConditionUnwrap）

让我们用聚类技术，自动定位6个集群。这里，使用的K-means算法。K-手段是最流行的聚类算法之一，但它可能会导致局部最优。这是一个很好的做法，重复K-均值聚类算法多次与不同的初始条件，并选择成本最低的结果。在这种情况下，该算法运行5次，结果是相同的。

选择= statset ('展示'，“最后一次”）;[clusterIndex, centers] = kmeans(table2array(opConditionUnwrap))， 6，．..“距离”，“sqeuclidean”，“复制”5,“选项”、选择);

重复1次迭代，总距离之和= 0.279547。重复2,1次迭代，总距离之和= 0.279547。重复3,1次迭代，总距离之和= 0.279547。重复4,1次迭代，总距离之和= 0.279547。重复5,1次迭代，总距离之和= 0.279547。距离的最佳总和= 0.279547

可视化聚类结果和识别的聚类中心。

图helperPlotClusters(opConditionUnwrap, clusterIndex, centers)

如图所示，聚类算法成功地找到了6种工作状态。

工作制度规范化

让我们对按不同工作机制分组的度量进行标准化。首先，计算每个传感器测量的平均值和标准偏差按上一节中确定的工作机制分组。

centerstats =结构('意思'，表（），“SD”，表（））；为v = dataVariables centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v))， clusterIndex);(@std, trainDataUnwrap.(char(v))) = splitapply(@std, trainDataUnwrap.(char(v))， clusterIndex);结束centerstats.Mean

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.92 1502.1 1311.4 10.52 15.493 394.32 23198784.1 1.26 45.496 371.44 2388.2 8133.9 8.6653 0.03 369.74 2319 100 28.525 17.115 518.67 642.71 1590.7 1409.4 14.62 21.61 553.3 2388.1 9062.3 1.3 47.557 521.36 2388.1 8140.9 8.4442 0.03 393.27 2388 100 38.809 23.285 462.54 536.87 1262.8 1050.6 7.05 9.0275 175.4 1915.4 8014.9 0.93989 36.808 164.56 2028.3 7878 10.916 0.02307.39 1915 84.93 14.262 8.5552 445 549.72 1354.7 1128.2 3.91 5.7158 138.62 2212 8327 1.0202 42.163 130.53 2388 8088.4 9.3775 0.02 331.12 2212 100 10.584 6.3515 491.19 607.59 1486 1253.6 9.35 13.657 334.46 2324 8729.1 1.0777 44.466 314.85 2388.2 8065 9.2347 0.022299 365.45 2324 100 24.447 14.67 449.44 555.82 1366.9 1131.9 5.48 8.0003 194.43 2223 8355.2 1.0203 41.995 183.01 2388.1 8071.1 9.3344 0.02 334.29 2223 100 14.827 8.8966

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  __________ ________ ________ ________ __________ _________ ________ ________ ________ __________ _________ _________ _________ _________ ___________________ _________ _________ __________ _________ _________ 1.4553e-11 0.47617 5.8555 8.3464 1.1618e-12 0.0047272 0.6536 0.094487 18.057 1.51e-13 0.25089 0.52838 0.096183 16.012 0.037598 9.9929e-16 1.4723 00 0.14522 0.086753 4.059e-11 0.48566 5.9258 8.8223 3.7307e-14 0.0011406 0.86236 0.068654 20.061 1.2103e-13 0.25937 0.71239 0.069276 17.21213.624 0.044142 8.6744e-16 1.3083 0 5.9833e-12 0.11212 0.06722 0 0.44169 1.9763e-13 0.0049401 0.44198 0.30713 18.389 0.0014951 0.23584 0.3432 0.33144 0.0371351.4174 00 0.10778 0.063991 4.4456e-11 0.46992 5.7664 7.8679 8.9892e-13 0.0046633 0.9984 0.13032 17.983 0.0042388 0.2362 0.9055 0.13285 15.792 0.038156 0.0042084 1.428 00 0.13352 0.0797310.069348

每个区域中的统计数据可用于规范化训练数据。对于每个集合成员，提取每行的操作点，计算其到每个簇中心的距离，并找到最近的簇中心。然后，对于每个传感器测量，减去平均值并除以该簇的标准偏差。如果偏差接近0时，将标准化传感器测量值设置为0，因为几乎恒定的传感器测量值对剩余使用寿命估计不有用。有关区域标准化功能的更多详细信息，请参阅最后一节“辅助功能”。

trainDataNormalized = cellfun(@(data)．..trainData,“UniformOutput”、假);

可视化根据工作状态规范化的数据。一些传感器测量的退化趋势在归一化后现现。

figure helperPlotEnsemble(trainDataNormalized, timeVariable, datavvariables (1:4)， nsample)

Trendability分析

现在，选择最具趋势的传感器测量值来构建一个用于预测的健康指标。对于每个传感器测量，估计线性退化模型，并对信号的斜率进行排序。

numSensors =长度(dataVariables);signalSlope = 0 (numSensors, 1);警告=警告(“关闭”）;为ct = 1:numSensors tmp = cellfun(@(tbl) tbl(:， cellstr(dataVariables(ct)))， traindatanorized，“UniformOutput”、假);mdl = linearDegradationModel ();%建立模型fit (mdl、tmp);%训练模式信号斜率（ct）=mdlθ；结束警告(警告);

排序的信号斜率并选择具有最大斜率8个传感器。

[~， idx] = sort(abs(signalSlope)，“下”）;sensorTrended =排序(idx (1:8))

sensorTrended =8×12 3 4 7 11 12 15 17

可视化选定的趋势传感器测量值。

figure helperPlotEnsemble(trainDataNormalized, timeVariable, datavvariables (sensorTrended(3:6))， nsample)

请注意，一些最trendable信号显示出积极的发展趋势，而其他人表现出消极趋势。

构建健康指示器

本节的重点熔合传感器测量到一个单一的健康指标，与其中基于相似性的模型被训练。

假设所有运行到故障的数据都以健康状态开始。开始时的健康状态分配值为1，故障时的健康状态分配值为0。健康状态假设随着时间的推移从1线性下降到0。此线性下降用于帮助融合传感器值。更多文献[2-5]中描述了复杂的传感器融合技术。

为j=1:numel(trainDataNormalized) data = trainDataNormalized{j};荷重软化= max (data.time) -data.time;数据。Health_condition = rul / max(rul);trainDataNormalized {j} =数据;结束

想象健康状况。

图helperPlotEnsemble (trainDataNormalized timeVariable,“健康状况”，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:结束)

重量=8×1-0.0308 -0.0308 -0.0536 0.0033 -0.0639 0.0051 -0.0408 -0.0382

通过将传感器测量值与相关权重相乘，构建单个运行状况指示器。

trainDataFused=cellfun（@（数据）降级传感器融合（数据、传感器使用、重量），trainDataNormalized，．..“UniformOutput”、假);

可视化训练数据的融合运行状况指示器。

图helperPlotEnsemble(trainDataFused， []， 1, nsample) xlabel(“时间”)伊拉贝尔(“健康指示器”)标题(“训练数据”）

来自多个传感器的数据融合到单个健康指示器中。健康指示器通过移动平均滤波器平滑。有关更多详细信息，请参阅上一节中的辅助函数“degradationSensorFusion”。

对验证数据应用相同的操作

对验证数据集重复状态归一化和传感器融合过程。

validationdatornormalize = cellfun(@(data))．..validationData,“UniformOutput”、假);validationDataFused = cellfun(@(data)降解传感器融合(data, sensorToFuse, weights)，．..validationDataNormalized,“UniformOutput”、假);

可视化验证数据的运行状况指示器。

图helperPlotEnsemble（validationDataFused，[]，1，NSample个）xlabel（“时间”)伊拉贝尔(“健康指示器”)标题(验证数据的）

建立相似规则模型

利用训练数据建立基于残差的相似度RUL模型。在这种情况下，模型试图用一个二阶多项式拟合每个融合数据。

数据间距离 $我$ 和数据 $j$ 由残差的1-范数计算

$d （我， j ）＝ | | y_{j} - {y_{}^{ˆ}}_{j ，我} | |_{1}$

哪里 $y_{j}$ 是机器的健康指标 $j$ ， ${y_{}^{ˆ}}_{j ，我}$ 估计的机器运行状况指示器是否正确 $j$ 使用机器所确定的2阶多项式模型 $我$ ．

相似性分数由以下公式计算

$年代 c o r e （我， j ）＝ e x p （ - d （我， j ）^{2} ）$

在验证数据集中给定一个集合成员，模型将在训练数据集中找到最接近的50个集合成员，并基于这50个集合成员拟合概率分布，并使用分布的中位数作为RUL的估计。

mdl = residualSimilarityModel (．..“方法”，“poly2”，．..“距离”，“绝对”，．..“NumNearestNeighbors”, 50岁,．..“标准化”1);fit (mdl trainDataFused);

绩效评估

为了评估相似度RUL模型，使用样本验证数据的50%、70%和90%来预测其RUL。

Breakpoint = [0.5, 0.7, 0.9];validationDataTmp = validationDataFused {3};%使用一个验证数据进行说明

在第一个断点之前使用验证数据，这是生存期的50%。

bpidx = 1;validationDataTmp50 = validationDataTmp(1:装天花板(结束*断点(bpidx)):);true = length(validationDataTmp) - length(validationDataTmp50);[estRUL, ciRUL, pdfRUL] = predictRUL(mdl, validationDataTmp50); / /确认数据

可视化截断为50%的验证数据及其最近的邻居。

图比较（mdl，validationDataTmp50）；

将估计的规则l与真实的规则l进行比较，并将估计的规则l的概率分布可视化。

图helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

有估计的RUL和真实RUL之间的相对大的误差当机器处于中间健康阶段。在这个例子中，最相似的10条曲线是靠近开头，但分叉当它们接近故障状态，导致在分配RUL大致两种模式。

在第二个断点之前使用验证数据，这是生命周期的70%。

bpidx = 2;validationDataTmp70 = validationDataTmp(1: cell (end*breakpoint(bpidx))，:);true = length(validationDataTmp) - length(validationDataTmp70);[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp70); / /确认数据图比较(mdl validationDataTmp70);

图helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

当观测到更多数据时，RUL估计得到增强。

在第三个断点(即生命周期的90%)之前使用验证数据。

bpidx = 3;validationDataTmp90 = validationDataTmp(1: cell (end*breakpoint(bpidx))，:);true = length(validationDataTmp) - length(validationDataTmp90);[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp90); / /确认数据图比较(mdl validationDataTmp90);

图helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

当机器接近故障时，RUL估计在本例中得到了更大的增强。

现在，对整个验证数据集重复相同的评估过程，计算每个断点的估计RUL和真实RUL之间的误差。

numValidation =长度(validationDataFused);numBreakpoint =长度(断点);错误= 0 (numValidation, numBreakpoint);为dataIdx = 1:numValidation tmpData = validationDataFused{dataIdx};为bpidx = 1:numBreakpoint tmpDataTest = tmpData(1: cell (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) hold从xlabel（“预测错误”)标题(“使用每个验证集合成员的前50%的RUL预测错误”)ax（2）=子批次（3,1,2）；持有在…上直方图(错误(:,2),“BinWidth”5,“归一化”，“pdf”) plot(x70, pdf70) hold从xlabel（“预测错误”)标题(“使用每个验证集合成员的前70%的RUL预测错误”) ax(3) = subplot(3,1,3);持有在…上直方图(错误(:,3),“BinWidth”5,“归一化”，“pdf”) plot(x90, pdf90) hold从xlabel（“预测错误”)标题(“使用每个验证集合成员的前90%的RUL预测错误”)连接轴（ax）

绘制的预测误差在一个盒子里的情节，以可视化的位数，位数25-75和异常值。

图箱线图（错误，“标签”,{“50%”，'70％”，“90%”}) ylabel (“预测错误”)标题('使用每个验证集合成员的不同百分比预测错误'）

计算和可视化预测误差的平均值和标准偏差。

errorMean =意味着(错误)

errorMean =1×3-5.8944 3.1359 3.3555

errorMedian =值(错误)

误差中值=1×3-4.8538 5.3763 3.6580

errorSD = STD（误差）

errorSD =1×326.4916 20.0720 18.0313

figure errorbar([50 70 90]， errorMean, errorSD，“o”，'MarkerEdgeColor'，“r”）XLIM（[40，100]）xlabel（“用于RUL预测的验证数据百分比”)伊拉贝尔(“预测错误”)传说(“具有1个标准差误差条的平均预测误差”，“位置”，“南”）

结果表明，随着观测数据的增加，误差更集中在0(更少的异常值)附近。

参考

[1] A. Saxena和K. Goebel(2008)。“PHM08挑战数据集”，NASA Ames预测数据存储库(http://ti.arc.nasa.gov/project/prognostic-data-repository)， NASA Ames研究中心，Moffett Field, CA

Roemer, Michael J.， Gregory J. Kacprzynski, Michael H. Schoeller。“利用健康管理信息融合改进诊断和预后评估。”AUTOTESTCON学报,2001。IEEE系统准备技术会议．IEEE 2001。

[3]格贝尔，楷，和皮耶罗Bonissone。“预后信息融合为恒定负载系统”。信息融合，2005年第八届信息技术国际会议. 第二卷。IEEE，2005年。

[4]王朋和David W.科伊特。“可靠性预测基于对多个降解措施系统退化建模。”可靠性和可维护性，2004年学术年会，RAMS．IEEE 2004。

[5] Jardine, Andrew KS, Daming Lin和Dragan Banjevic。“对机械诊断和基于状态维护的预测进行了回顾。”机械系统和信号处理20.7(2006): 1483 - 1510。

辅助函数

函数规范化(data, centers, centerstats)%对数据的每个观察(行)执行归一化%根据观察所属的群集。conditionIdx = 3：5;dataIdx = 6:26;%行操作data{:， dataIdx} = table2array(．..rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats)，．..数据,“SeparateInputs”、假));结束函数rownormalization = localNormalize(row, conditionIdx, dataIdx, centers, centerstats)%每行的规范化操作。％获取工作点和传感器测量OPS =行（1，conditionIdx）;传感器=行（1，dataIdx）;％查找哪个集群中心是最接近Dist = sum((centers - ops)。2) ^ 2;[~， idx] = min(dist);％归一化的传感器测量由群集的平均值和标准偏差。将NaN和Inf重新赋值为0。rowNormalized = (sensor - centerstats。Mean{idx，:}) ./ centerstats。SD {idx:};isnan(rowNormalized) | isinf(rowNormalized)) = 0;结束函数dataFused =降解sensorfusion (data, sensorToFuse, weights)％组合来自不同，根据传感器的测量%的权重，平滑融合数据并偏移数据%，以便所有数据从1开始%根据权重融合数据dataToFuse = data{:， cellstr(sensorToFuse)};dataFusedRaw = dataToFuse *权重;％光滑移动平均值的融合的数据stepBackward = 10;stepForward = 10;dataFused = movmean(dataFusedRaw， [stepback stepForward]); / /将数据还原%将数据偏移到1dataFused = dataFused + 1 - dataFused(1);结束