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

打开直播脚本

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

数据准备

由于数据集很小，因此可以将整个降级数据加载到内存中。从中下载并解压缩数据集https://ti.art.nasa.gov/tech/dash/groups/pcoe/prognostic-data-repository/到当前目录。使用helperLoadData辅助功能要加载并将培训文本文件转换为时间表的单元格数组。训练数据包含218个运行失败模拟。这组测量被称为“合奏”。

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

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

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

头(degradationData {1})

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 __ ____ ____________ ____________ ____________ ________ ________ ________ ________ ________ ________ ________ ________ ________ _________ _________ _________ _________ _________ _________ _________ _________ _________ __________________ _________ 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 3912388 100 38.99 9 0.8401343.9222.03085/303855/30123557 0.706/303858034.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.99.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.081589.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 (“默认”）%以确保结果的可重复性numemble =长度（degraationdata）;numfold = 5;cv = cvpartition（数值，“KFold”, numFold);trainData =降解数据(training(cv, 1));validationData =降解数据(测试(cv, 1));

指定感兴趣的变量组。

varnames = string（degadationdata {1} .properties.variablenames）;timevariable = varnames {2};条件variables = varnames（3：5）;datavariables = varnames（6:26）;

可视化集合数据的样本。

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（TrainData {：}）;opconditionUnwrap = traindataUnwrap（：，cellstr（CitnitionVariables））;

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

图辅助全架子（OpconditionUnwrap）

让我们使用集群技术来自动定位这6个集群。这里使用K-means算法。K-means是最流行的聚类算法之一，但它会导致局部最优。在不同的初始条件下多次重复K-means聚类算法，选取代价最小的结果是一个很好的实践。在这种情况下，算法运行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

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

Figure Helperplotclusters（OpconditionUnwrap，ClusterIndex，中心）

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

工作制度正常化

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

Centerstats = struct（“的意思是”， 桌子（），“SD”， 桌子（））;为v = dataVariables centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v))， clusterIndex);(@std, trainDataUnwrap.(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.92 1502.1 1311.4 10.52 15.493 394.32 2319 8784.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.81050.6 7.05 9.0275 175.4 1915.4 8014.9 0.93989 36.808 164.56 2028.3 7878 445 1915 10.916 0.02 307.39 84.93 14.262 8.5552 549.72 1354.7 1128.2 3.91 5.7158 138.62 2212 8327 2388 1.0202 42.163 130.53 8088.4 9.3775 0.02 331.12 2212 1486 100 10.584 6.3515 491.19 607.59 1253.6 9.35 13.657 334.46 2324 8729.1 1.0777 44.466 314.85 2388.2 8065 9.23470.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.65536 0.094487 18.057 1.51e-13 0.25089 0.52838 0.096183 16.012 0.037598 9.9929e-16 1.4723 0 0 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.212 0.036504 9.9929E-16 1.5046 0 0 0.17565 0.10515 2.7685e-11 0.35468 5.2678 6.9664 2.043E-14 0.0043301 0.45074 0.2743 14.741 0.0010469 0.21539 0.33985 0.28932 13.624 0.044142 8.6744e-16 1.3083 0 5.9833e-12 0.11212 0.06722 0 0.44169 5.6853 7.5741 1.9763e-13 0.00494010.44198 0.30713 18.389 0.0014951 0.23584 0.3432 0.33144 16.917 0.037135 2.231e-15 1.4174 0 0 0.10778 0.063991 4.4456e-11 0.46992 5.7664 7.8679 8.9892e-13 0.0046633 0.59984 0.13032 17.983 0.0042388 0.2362 0.49055 0.13285 15.792 0.038156 0.0042084 1.4428 0 0 0.13352 0.079731 4.3489e-11 0.44341 5.7224 7.4842 3.7485e-13 0.0017642 0.4734 0.28889 17.608 0.0017978 0.23242 0.38243 0.31 15.91 0.038624 8.5009E-16 1.3974 0 0 0.11311 0.069348

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

trainDataNormalized = cellfun(@(data)．..trainData,'统一输出'、假);

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

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

Trendability分析

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

numsensors =长度（datavariables）;信号LOPE =零（NUMSENSORS，1）;warn =警告（“关闭”）;为ct = 1:numSensors tmp = cellfun(@(tbl) tbl(:， cellstr(dataVariables(ct)))， traindatanorized，'统一输出'、假);mdl = linearDegradationModel ();%建立模型适合（MDL，TMP）;%训练模式信号升（CT）= mdl.theta;结束警告（警告）;

对信号斜率进行排序，选择斜率最大的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)

请注意，一些最具趋势性的信号显示出积极的趋势，而另一些则显示出消极的趋势。

构建健康指示器

本节重点介绍将传感器测量结果融合到单个健康指标中，使用该指标训练基于相似性的模型。

假设所有运行失败的数据都以健康的状态开头。开始时的健康状况为1的值1，并且失败时的健康状况被分配了值0.假设健康状况随时间从1到0线性降级。这种线性劣化用于帮助熔断传感器值。文献中描述了更复杂的传感器融合技术[2-5]。

为j = 1：numel（traindatanormalized）data = traindatanormalized {j};rul = max（data.time）-data.time;data.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.coffients.imate（2：结束）

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

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

TrainDatafused = Cellfun（@（数据）降级病症（数据，sensortofuse，重量），traindatanormization，．..'统一输出'、假);

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

Figure Helperplotysemble（TrainDatafused，[]，1，Nsample）Xlabel（“时间”）ylabel（“健康指示器”)标题(的训练数据）

来自多个传感器的数据融合到单个健康指示器中。健康指示器由移动平均滤波器平滑。有关更多详细信息，请参阅最后一节中的帮助函数“降级损坏”。

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

使用验证数据集重复调节归一化和传感器融合过程。

validationdatornormalize = cellfun(@(data))．..validationData,'统一输出'、假);ValidationDatafused = CellFun（@（数据）降级损坏（数据，Sensortofuse，权重），．..验证atanormalized，'统一输出'、假);

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

图helperPlotEnsemble(validationDataFused， []， 1, nsample)“时间”）ylabel（“健康指示器”)标题(验证数据的）

建立相似规则模型

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

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

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

在哪里 $y_{j}$ 机器的运行状况指示器是否正常 $j$ ， ${y_{}^{ˆ}}_{j ，我}$ 估计的机器运行状况指示器是否正确 $j$ 利用机器辨识的二阶多项式模型 $我$ ．

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

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

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

MDL =残留匹配性模型（．..'方法'，'poly2'，．..'距离'，'绝对'，．..'numnearestneighbors', 50岁,．..“标准化”1);fit (mdl trainDataFused);

绩效评估

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

断点= [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] = ksdites（错误（：，1））;[pdf70，x70] = ksdites（错误（:,2））;[pdf90，x90] = ksdites（错误（:,3））;图AX（1）=子图（3,1,1）;持有在直方图(错误(:1),“BinWidth”5,“归一化”，“pdf”）绘图（x50，pdf50）持有离开包含(预测误差的)标题(“使用每个验证集合成员的前50%的RUL预测错误”）斧头（2）=子图（3,1,2）;持有在直方图(错误(:,2),“BinWidth”5,“归一化”，“pdf”) plot(x70, pdf70) hold离开包含(预测误差的)标题('RUL预测错误使用每个验证集合成员的前70％') ax(3) = subplot(3,1,3);持有在直方图(错误(:,3),“BinWidth”5,“归一化”，“pdf”) plot(x90, pdf90) hold离开包含(预测误差的)标题('RUL预测错误使用每个验证集合成员的前90％'）Linkaxes（AX）

在框图中绘制预测误差，以可视化中值、25-75分位数和离群值。

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

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

errorMean =意味着(错误)

errormean =1×3.-5.8944 3.1359 3.3555

errormedian =中位数（错误）

errormedian =1×3.-4.8538 5.3763 3.6580

errorSD =性病(错误)

errorSD =1×3.26.4916 20.0720 18.0313

figure errorbar([50 70 90]， errorMean, errorSD，“o”，“MarkerEdgeColor”，“r”) xlim([40,100])“用于RUL预测的验证数据百分比”）ylabel（预测误差的)传说(“具有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。

戈贝尔，凯，皮耶罗·博尼松。“恒定负荷系统的预测信息融合。”信息融合，2005年第8届国际会议．卷。2. IEEE，2005。

Wang, Peng和David W. Coit。“基于退化建模的具有多种退化措施的系统可靠性预测。”可靠性和可维护性，2004年年度研讨会．IEEE，2004年。

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

辅助函数

功能数据= Regimenormalization（数据，中心，Centerstats）%对数据的每个观察(行)执行归一化%根据观察所属的群集。conditionIdx = 3:5;dataIdx =第一;%行操作数据{：，dataidx} = table2array（．..rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats)，．..数据,“SeparateInputs”、假));结束功能rownormalization = localNormalize(row, conditionIdx, dataIdx, centers, centerstats)每行的％归一化操作。获得操作点和传感器测量值ops = row(1, conditionIdx);sensor = row(1, dataIdx); / /读取数据找出哪个集群中心是最近的dist = sum（（中心 -  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 *权重;%用移动平均值平滑融合数据STAPBACKWARD = 10;stepforward = 10;datafused = movmean（datafusedraw，[stepbackward stepforward]）;％偏移数据到1datafused = datafused + 1  -  datafused（1）;结束