该示例显示了如何构建指数劣化模型,以实时预测风力涡轮机轴承的剩余使用寿命(RUL)。指数劣化模型基于其参数前导者和最新测量预测RUL(历史记录到故障数据可以帮助估计模型参数前沿,但不需要它们)。该模型能够实时检测显着的降级趋势,并在新的观察变得可用时更新其参数前沿。该示例介绍了典型的预测工作流程:数据导入和探索,功能提取和后处理,特征重要性排名和融合,模型拟合和预测,以及性能分析。据/p>
数据集由由20齿小齿轮齿轮驱动的2MW风力涡轮机高速轴收集[1]。每天在连续50天中获得6秒的振动信号(3月17日有2个测量,在该示例中被视为两天)。内部竞争故障开发并导致轴承失效,在50天期间。据/p>
工具箱中提供了数据集的压缩版本。要使用压缩数据集,请将数据集复制到当前文件夹并启用其写入权限。据/p>
拷贝文件(据span style="color:#0000FF">......据/span>fullfile(matlabroot,据span style="color:#A020F0">'工具箱'据/span>那据span style="color:#A020F0">'predmaint'据/span>那据span style="color:#0000FF">......据/span>'predmaintdemos'据/span>那据span style="color:#A020F0">'windTurbineHighSpeedBearingPrognosis'据/span>),据span style="color:#0000FF">......据/span>'WindTurbineHighSpeedBearingPrognosis-DATA-硕士据/span>fileattrib(fullfile(据span style="color:#A020F0">'WindTurbineHighSpeedBearingPrognosis-DATA-硕士据/span>那据span style="color:#A020F0">'*。垫'据/span>),据span style="color:#A020F0">'+ w'据/span>)据/pre>
Compact DataSet的测量时间步长为5天。据/p>
TimeUnit =.据span style="color:#A020F0">'\ times 5天'据/span>;据/pre>
对于整个数据集,请访问以下链接据a href="https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data" target="_blank">https://github.com/mathworks/windturbinehighspeedbearingprognosiss-data.据/a>,将整个存储库作为zip文件下载,并将其保存在与此实时脚本相同的目录中。使用此命令解压缩文件。完整数据集的测量时间步长为1天。据/p>
如果据/span>存在(据span style="color:#A020F0">'windturbinehighpeedbearingprognosis-data-master.zip'据/span>那据span style="color:#A020F0">'文件'据/span>)解压缩(据span style="color:#A020F0">'windturbinehighpeedbearingprognosis-data-master.zip'据/span>)TIMEUNIT =据span style="color:#A020F0">'日'据/span>;据span style="color:#0000FF">结尾据/span>
在这个实施例的结果是从完整数据集生成的。强烈建议下载完整数据集来运行这个例子。从小型数据集生成的结果可能没有意义。据/p>
创建一个据code class="literal">fillensembledataStore.据/code>的风力涡轮机数据。数据中包含的振动信号和转速计信号。这据code class="literal">fillensembledataStore.据/code>将解析文件名并将日期信息提取为据code class="literal">IndependentVariables据/code>。有关更多详细信息,请参阅与此示例相关联的支持文件中的辅助功能。万博1manbetx据/p>
hsbearing = fileEnsembleDatastore(据span style="color:#0000FF">......据/span>完整文件(据span style="color:#A020F0">'。'据/span>那据span style="color:#A020F0">'WindTurbineHighSpeedBearingPrognosis-DATA-硕士据/span>),据span style="color:#0000FF">......据/span>'。垫'据/span>);hsbearing.datavariables = [据span style="color:#A020F0">“振动”据/span>那据span style="color:#A020F0">“tach”据/span>]; hs.独立变量=据span style="color:#A020F0">“日期”据/span>; hsbearing.SelectedVariables=[据span style="color:#A020F0">“日期”据/span>那据span style="color:#A020F0">“振动”据/span>那据span style="color:#A020F0">“tach”据/span>];hsbearing.readfcn = @helperreaddata;hsbearing.writetomemberfcn = @helperwriteTohsbearing;高大(Hsbearing)据/pre>
ANS = M×3高大表日期振动转速____________________ _________________ _______________ 07-MAR-2013 1时57分46秒[585936×1双] [2446×1双] 08-MAR-2013 2点34分21秒[585936×1双] [2411×1双] 09-MAR-2013二时33分43秒[585936×1双] [2465×1双] 10-MAR-2013 3点01分02秒[585936×1双] [2461×1双] 11-MAR-2013 3时00分24秒[585936×1双] [2506×1双] 12-MAR-2013六时17分10秒[585936×1双] [2447×1双] 13-MAR-20136点34分04秒[585936×1双] [2438×1双] 14-MAR-2013 6点五十分41秒[585936×1双] [2390×1双]::::::据/pre>
振动信号的采样率为97656 Hz。据/p>
FS = 97656;据span style="color:#228B22">%Hz.据/span>
本节探讨了时域和频域中的数据,并寻求启发用于提取预后目的的功能。据/p>
首先在时域中可视化振动信号。在该数据集中,在连续50天中有50个振动信号6秒。现在绘制50振动信号彼此一个。据/p>
重置(hsbearing)tstart = 0;图持有据span style="color:#A020F0">在据/span>尽管据/span>hasdata(hsbearing)数据=读取(hsbearing);v = data.vibration {1};t = tstart +(1:长度(v))/ fs;据span style="color:#228B22">%下采样信号以减少存储器使用据/span>绘图(t(1:10:end),v(1:10:end))t开始=t(end);据span style="color:#0000FF">结尾据/span>抓住据span style="color:#A020F0">离开据/span>Xlabel(据span style="color:#A020F0">'时间,每天6秒,总共50天'据/span>)ylabel(据span style="color:#A020F0">'加速(g)'据/span>)据/pre>
在时域上的振动信号揭示了信号冲动的增加趋势。指标量化信号的冲动,如峰度,峰 - 峰值,波峰因数据span class="emphasis">等等。据/em>,对于该风力涡轮机轴承数据集进行潜在的预后特征[2]。据/p>
在另一方面,光谱峰度被认为是在频域[3]风力涡轮机的预后有力工具。为了显现沿着时间的光谱峰度的变化,绘制光谱峰度值作为频率的函数和测量一天。据/p>
颜色栏中指示的故障严重性是测量日期为0到1比例。观察到,随着机器条件降解,速度约为10kHz的光谱峰值值逐渐增加。光谱峰度的统计特征,如平均值,标准偏差据span class="emphasis">等等。据/em>,将是轴承劣化的潜在指标[3]。据/p>
基于前一节中的分析,将提取源自时域信号和光谱峰度的统计特征的集合。[2-3]中提供了有关该功能的更多数学细节。据/p>
首先,在将其写入FileSeneMbleyataStore之前,在DataVaria中预先分配特征名称。据/p>
对于每个集合构件计算特征值。据/p>
选择自变量据code class="literal">日期据/code>并将提取的所有特征构建特征表。据/p>
由于要素表足够小以适合内存(50到15),因此在处理前收集表。对于大数据,建议以高表格执行操作,直到您相信输出足够小以适合内存。据/p>
将表转换为时间表,以便时间信息始终与特征值相关联。据/p>
提取的特征通常与噪声相关联。具有相反趋势的噪音有时可能对rul预测有害。另外,要引入的特征性能指标,单调性待引入的一个是不稳健的噪声。因此,将具有5个步骤的滞后窗口的因果移动平均滤波器应用于提取的特征,其中“因果”表示移动平均滤波中没有将未来的值用于未来的值。据/p>
这是一个示例,示出了平滑之前和之后的功能。据/p>
移动平均平滑引入了信号的时间延迟,但是通过在RUL预测中选择适当的阈值,可以减轻延迟效果。据/p>
在实践中,在开发预后算法时,整个生命周期的数据不可用,但是假设已经收集了生命周期早期的一些数据是合理的。因此,在前20天内收集的数据(40%的生命周期)被视为培训数据。以下特征重要性排名和融合仅基于培训数据。据/p>
在该示例中,[3]提出的单调性用于量化预后目的的特征的优点。据/p>
的单调性据span class="inlineequation">
特征据span class="inlineequation">
计算为据/p>
在哪里据span class="inlineequation">
是测量点数,在这种情况下据span class="inlineequation">
。据span class="inlineequation">
在机器的数量监测,在这种情况下,据span class="inlineequation">
。据span class="inlineequation">
是个据span class="inlineequation">
测量的功能据span class="inlineequation">
机器。据span class="inlineequation">
,即信号的差异据span class="inlineequation">
。据/p>
该信号的峰度是基于单调顶部的特征。据/p>
选择具有大于0.3的特征重要性分数的功能对于下一节中的功能融合。据/p>
主成分分析(PCA)用于该示例中的尺寸减小和特征融合。在执行PCA之前,它是一个很好的做法,以将特征标准化为相同的规模。注意,PCA系数和归一化中使用的平均值和标准偏差是从训练数据获得的,并应用于整个数据集。据/p>
平均值,标准偏差和PCA系数用于处理整个数据集。据/p>
可视化前两个主组件的空间中的数据。据/p>
该图表明,随着机器接近失败第一主成分正在增加。因此,第一主成分是一种很有前途稠健康指示符。据/p>
可视化健康指标。据/p>
指数劣化模型被定义为据/p>
在哪里据span class="inlineequation">
是健康指标作为时间的函数。据span class="inlineequation">
被认为是一个常数截距项。据span class="inlineequation">
和据span class="inlineequation">
是决定模型斜率的随机参数,其中据span class="inlineequation">
是Lognormal分布的据span class="inlineequation">
是高斯分布的。在每一步据span class="inlineequation">
,分布据span class="inlineequation">
和据span class="inlineequation">
更新为基础的最新观测后据span class="inlineequation">
。据span class="inlineequation">
是高斯白噪声屈服于据span class="inlineequation">
. 这个据span class="inlineequation">
指数中的术语是预期据span class="inlineequation">
满足据/p>
。据/p>
这里,指数劣化模型适合最后一节中提取的健康指示符,并且在下一节中评估性能。据/p>
首先移动运行状况指示器,使其从0开始。据/p>
阈值的选择通常基于机器的历史记录或某种域特定知识。由于此数据集中未使用历史数据,因此选择健康指示符的最后值作为阈值。建议基于平滑(历史)数据选择阈值,以便部分减轻平滑的延迟效果。据/p>
如果历史数据是可用的,使用据code class="literal">合身据/code>提供的方法据code class="literal">exponentialDegradationModel据/code>估计先验和截距。然而,历史数据是无法用于此风力涡轮机轴承的数据集。现有斜率参数的被选择任意大方差(据span class="inlineequation">
)因此该模型主要依赖于观察到的数据。基于据span class="inlineequation">
,拦截据span class="inlineequation">
被设定为据span class="inlineequation">
因此,模型也将从0开始。据/p>
健康指标变异与噪声变异之间的关系可以推导为据/p>
这里假设噪声的标准偏差在接近阈值时导致健康指示器的10%变化。因此,噪声的标准偏差可以表示为据span class="inlineequation">
。据/p>
指数衰减模型还提供了一个功能,以评估斜坡的意义。一旦检测到健康指示器的显著斜率,该模型会忘记以前的意见,然后重新启动在原有基础上先验估计。检测算法的灵敏度可以通过指定被调谐据code class="literal">SlopeDetectionLevel据/code>。如果P值小于据code class="literal">SlopeDetectionLevel据/code>,声明斜率被检测到。这里据code class="literal">SlopeDetectionLevel据/code>设置为0.05。据/p>
现在使用上面讨论的参数创建指数劣化模型。据/p>
用据code class="literal">预测据/code>和据code class="literal">更新据/code>方法预测RUL并实时更新参数分布。据/p>
这里是实时RUL估计的动画。据/p>
-据span class="inlineequation">
情节用于预后性能分析[5],其中据span class="inlineequation">
绑定设定为20%。估计的rul之间的概率在于据span class="inlineequation">
真正的RUL的约束计算为模型的性能指标:据/p>
在哪里据span class="inlineequation">
是估计的统治据span class="inlineequation">
那据span class="inlineequation">
是真正的rul据span class="inlineequation">
那据span class="inlineequation">
是估计的模型参数吗?据span class="inlineequation">
。据/p>
由于预设事先不反映真实的先前,模型通常需要几个时间步骤来调整到适当的参数分布。随着更多数据点可用,预测变得更加准确。据/p>
可视化预测rul内的概率据span class="inlineequation">
跳跃据/p>
[1] Bechhoefer,Eric,Brandon Van Hecke和David他。“改善光谱分析的处理。”据span class="emphasis">预测和健康管理协会,新奥尔良,洛杉矶,洛杉矶的年会据/em>。2013年。据/p>
[2] Ali,Jaouher Ben,等。“在无监督机器学习的实验条件下,在线自动诊断风力涡轮机轴承逐步降解。”据span class="emphasis">应用声学据/em>132 (2018): 167-181.据/p>
[3]赛迪,卢特菲,等。“风力涡轮机高速通过的光谱峰度衍生指数和SVR轴轴承健康预后。”据span class="emphasis">应用声学据/em>120(2017):1-8。据/p>
[4]积分,杰米baalis。“合并数据来源以预测剩余的使用寿命 - 一种识别预后参数的自动化方法。”(2010)。据/p>
[5] Saxena先生,阿比纳夫,等。“指标的预后性能下线的评价。”据span class="emphasis">国际杂志预测与健康管理据/em>1.1(2010):4-23。据/p>
hsbearing.datavariables = [据span style="color:#A020F0">“振动”据/span>那据span style="color:#A020F0">“tach”据/span>那据span style="color:#A020F0">“SpectralKurtosis”据/span>];颜色= parula(50);图持有据span style="color:#A020F0">在据/span>复位(hsbearing)天= 1;据span style="color:#0000FF">尽管据/span>hasdata(hsbearing)数据=读取(hsbearing);data2add =表;据span style="color:#228B22">%获得振动信号和测量日期据/span>v = data.vibration {1};据span style="color:#228B22">%计算谱峰度与窗口大小= 128据/span>WC = 128;[SK,F] = pkurtosis(V,FS,WC);data2add.SpectralKurtosis = {表(F,SK)};据span style="color:#228B22">%绘制光谱峰度据/span>Plot3(F,Day * Ones(尺寸(f)),sk,据span style="color:#A020F0">'颜色'据/span>,颜色(日,:))据span style="color:#228B22">%写光谱峰值值据/span>WriteTolastmemberRead(Hsbearing,Data2Add);据span style="color:#228B22">%递增天数据/span>天=天+ 1;据span style="color:#0000FF">结尾据/span>抓住据span style="color:#A020F0">离开据/span>Xlabel(据span style="color:#A020F0">'频率(Hz)'据/span>)ylabel(据span style="color:#A020F0">“时间(天)”据/span>)Zlabel(据span style="color:#A020F0">'光谱峰氏'据/span>) 网格据span style="color:#A020F0">在据/span>查看(-45,30)cbar =彩色杆;ylabel(c bar,据span style="color:#A020F0">'故障严重性(0 - 健康的,1 - 故障)'据/span>)据/pre>
特征提取据/h3>
hsbearing.datavariables = [hsbearing.datavariables;据span style="color:#0000FF">......据/span>“意思”据/span>;据span style="color:#A020F0">“std”据/span>;据span style="color:#A020F0">“歪斜”据/span>;据span style="color:#A020F0">“kurtosis”据/span>;据span style="color:#A020F0">“Peak2Peak”据/span>;据span style="color:#0000FF">......据/span>“rms”据/span>;据span style="color:#A020F0">“crestfactor”据/span>;据span style="color:#A020F0">“ShapeFactor”据/span>;据span style="color:#A020F0">“暂停”据/span>;据span style="color:#A020F0">“marginfactor”据/span>;据span style="color:#A020F0">“活力”据/span>;据span style="color:#0000FF">......据/span>“skmean”据/span>;据span style="color:#A020F0">“skstd”据/span>;据span style="color:#A020F0">“倾斜度”据/span>;据span style="color:#A020F0">“skkurtosis”据/span>];据/pre>
hsbearing.SelectedVariables = [据span style="color:#A020F0">“振动”据/span>那据span style="color:#A020F0">“SpectralKurtosis”据/span>];重置(Hsbearing)据span style="color:#0000FF">尽管据/span>hasdata(hsbearing)数据=读取(hsbearing);v = data.vibration {1};sk = data.spectralkurtosis {1} .sk;特点=表;据span style="color:#228B22">%时域特征据/span>特征.. anean =均值(v);特征.std = std(v);特征.skewness = skewness(v);特征。kurtosis = kurtosis(v);特征.peak2peak = peak2peak(v);特征.rms = rms(v);特征.CrestFactor = Max(v)/ features.rms;特征.Shapefactor =特征.rms /均值(abs(v));特征.Impulsefactor = max(v)/平均值(abs(v));特征.marginfactor = max(v)/平均值(abs(v))^ 2; features.Energy = sum(v.^2);据span style="color:#228B22">%光谱峰峰相关特征据/span>特征.skmean =平均(sk);特征.skstd = std(sk);特征.skskewness = skewness(sk);特征.skkurtosis = kurtosis(sk);据span style="color:#228B22">%写导出的特征为相应的文件据/span>WriteTolastMemberRead(Hsbearing,功能);据span style="color:#0000FF">结尾据/span>
hsbearing.SelectedVariables = [据span style="color:#A020F0">“日期”据/span>那据span style="color:#A020F0">“意思”据/span>那据span style="color:#A020F0">“std”据/span>那据span style="color:#A020F0">“歪斜”据/span>那据span style="color:#A020F0">“kurtosis”据/span>那据span style="color:#A020F0">“Peak2Peak”据/span>那据span style="color:#0000FF">......据/span>“rms”据/span>那据span style="color:#A020F0">“crestfactor”据/span>那据span style="color:#A020F0">“ShapeFactor”据/span>那据span style="color:#A020F0">“暂停”据/span>那据span style="color:#A020F0">“marginfactor”据/span>那据span style="color:#A020F0">“活力”据/span>那据span style="color:#0000FF">......据/span>“skmean”据/span>那据span style="color:#A020F0">“skstd”据/span>那据span style="color:#A020F0">“倾斜度”据/span>那据span style="color:#A020F0">“skkurtosis”据/span>];据/pre>
featureTable =聚集(高(hsbearing));据/pre>
使用并行池“本地”评估高表达: - PASS 1为1:1秒评估完成1秒据/pre>
featureTable = table2timetable(featureTable)据/pre>
featureTable =据span class="emphasis">50×15时间表据/em>日均值标准偏度峰度Peak2Peak RMS峰值因子形状因子ImpulseFactor MarginFactor能源SKMean SKStd SKSkewness SKKurtosis ____________________ _______ ______ ___________ ________ _________ ______ ___________ ___________ _____________ ____________ __________ __________ ________ __________ __________ 07-MAR-2013 1时57分46秒0.34605 2.2705 0.0038699 2.9956 21.621 2.2967 4.9147 1.25356.1607 3.3625 3.0907e + 06 0.001253 0.025674 -0.22882 3.362 08-MAR-2013 2时34分21秒0.24409 2.0621 0.0030103 3.0195 19.31 2.0765 4.9129 1.2545 6.163 3.7231 2.5266e + 06 0.0046823 0.020888 0.057651 3.3508 09-MAR-2013 2点33分43秒0.218732.1036 -0.0010289 3.0224 21.474 2.1149 5.2143 1.2539 6.5384 3.8766 2.6208e + 06 -0.0011084 0.022705 -0.50004 4.9953 10-MAR-2013 3点01分02秒0.21372 2.0081 0.001477 3.0415 19.52 2.0194 5.286 1.2556 6.637 4.1266 2.3894e + 06 0.0087035 0.034456 1.4705 8.1235 11-MAR-2013三点00分24秒0.21518 2.0606 0.0010116 3.0445 21.217 2.0718 5.0058 1.2554 6.2841 3.8078 2.515e + 06 0.013559 0.032193 0.11355 3.848 12-MAR-2013 6点17分10秒0.29335 2.0791 -0.008428 3.018 20.05 2.0997 4.7966 1.2537 6.0136 3.5907 2.5833e + 06 0.0017539 0.028326 -0.1288 3.8072 13-MAR-2013 6点34分04秒0.21293 1.972 -0.0014294 3.0174 18.837 1.9834 4.84961.2539 6.0808 3.8441 2.3051e + 06 0.0039353 0.029292 -1.4734 8.1242 14-MAR-2013 6点五十分41秒0.24401 1.8114 0.0022161 3.0057 17.862 1.8278 4.8638 1.2536 6.0975 4.1821 1.9575e + 06 0.0013107 0.022468 0.27438 2.8597 15-MAR-2013 6时五十分03秒0.20984 1.9973 0.001559 3.0711 21.12 2.0083 5.4323 1.2568 6.8272 4.2724 2.3633e + 06 0.0023769 0.047767 -2.5832 20.171 16-MAR-2013六点56分43秒0.23318 1.9842 -0.0019594 3.0072 18.832 1.9979 5.0483 1.254 6.3304 3.9734 2.3387e + 06 0.0076327 0.030418 0.52322 4.0082 17-MAR-2013 6点56分04秒0.21657 2.113 -0.0013711 3.1247 21.858 2.1241 5.4857 1.2587 6.9048 4.0916 2.6437e + 06 0.0084907 0.037876 2.3753 11.491 17-MAR-2013十八时47分56秒0.19381 2.1335 -0.012744 3.0934 21.589 2.1423 4.7574 1.2575 5.9825 3.5117 2.6891e+06 0.019624 0.055537 3.1986 17.796 18-MAR-2013 18:47:15 0.21919 2.1284 -0.0002039 3.1647 24.051 2.1396 5.7883 1.2595 7.2902 4.2914 2.6824e + 06 0.016315 0.064516 2.8735 11.632 20-MAR-2013〇时33分54秒0.35865 2.2536 -0.002308 3.0817 22.633 2.2819 5.2771 1.2571 6.6339 3.6546 3.0511e + 06 0.0011097 0.051539 -0.0567747.0712 21-MAR-2013零时33分14秒0.1908 2.1782 -0.00019286 3.1548 25.515 2.1865 6.0521 1.26 7.6258 4.3945 2.8013e + 06 0.0040392 0.066254 -0.39587 12.111 22-MAR-2013 0时39分50秒0.20569 2.1861 0.0020375 3.2691 26.439 2.1958 6.1753 1.2633 7.80114.4882 2.825e + 06 0.020676 0.077728 2.6038 11.088⋮据/pre>
功能后处理据/h3>
variableNames = featureTable.Properties.VariableNames;featureTableSmooth = varfun(@(X)movmean(X,[5 0]),featureTable);featureTableSmooth.Properties.VariableNames = variableNames;据/pre>
图持有据span style="color:#A020F0">在据/span>图(featureTable.Date,featureTable.SKMean)图(featureTableSmooth.Date,featureTableSmooth.SKMean)保持据span style="color:#A020F0">离开据/span>Xlabel(据span style="color:#A020F0">'时间'据/span>)ylabel(据span style="color:#A020F0">'特征值'据/span>) 传奇(据span style="color:#A020F0">“以前平滑”据/span>那据span style="color:#A020F0">'平滑后'据/span>) 标题(据span style="color:#A020F0">'skmean'据/span>)据/pre>
培训数据据/h3>
breaktime = datetime(2013,3,27);breakpoint = find(featureTablesmooth.date
功能重要性排序据/h3>
由于移动窗口平滑已经完成,因此将“Windowsize”设置为0到0据/span>%关闭功能内的平滑据/span>FeatureImportance =单调性(TrainData,据span style="color:#A020F0">'Windowsize'据/span>,0);helperSortedBarPlot(featureImportance,据span style="color:#A020F0">“单调”据/span>);据/pre>
TrainDataselected = TrainData(:, FeedingImportance {:,:}> 0.3);feationElted = featureTablesmooth(:, FeedingImportance {:,:}> 0.3)据/pre>
feedingeled =据span class="emphasis">50×5时间表据/em>日均值峭度形状因子MarginFactor SKStd ____________________ _______ ________ ___________ ____________ ________ 07-MAR-2013 1时57分46秒0.34605 2.9956 1.2535 3.3625 0.025674 08-MAR-2013 2时34分21秒0.29507 3.0075 1.254 3.5428 0.023281 09-MAR-2013 02:33:43 0.26962 3.0125 1.254 3.6541 0.023089 10-MAR-2013 3点01分02秒0.25565 3.0197 1.2544 3.7722 0.025931 11-MAR-2013 3点00分24秒0.24756 3.0247 1.2546 3.7793 0.027183 12-MAR-2013六时17分10秒0.25519 3.0236 1.25443.7479 0.027374 13-MAR-2013 6点34分04秒0.233 3.0272 1.2545 3.8282 0.027977 14-MAR-2013六时50分41秒0.23299 3.0249 1.2544 3.9047 0.02824 15-MAR-2013 6点50分03秒0.2315 3.033 1.2548 3.9706 0.032417 3月16日-2013 6点56分43秒0.23475 3.0273 1.2546 3.9451 0.031744 17-MAR-2013 6点56分04秒0.23498 3.0407 1.2551 3.9924 0.032691 17-MAR-2013 18时47分56秒0.21839 3.0533 1.2557 3.9792 0.037226 18-MAR-2013 18点47分:15 0.21943 3.0778 1.2567 4.0538 0.043097 20-MAR-2013○点33分54秒0.23854 3.0905 1.2573 3.9658 0.047942 21-MAR-2013○点33分14秒0.23537 3.1044 1.2578 3.9862 0.05102322-MAR-2013零时39分50秒0.23079 3.1481 1.2593 4.072 0.058908⋮据/pre>
尺寸减少和特征融合据/h3>
含义=卑鄙(TrainDataselected {:,:,:});sdtrain = std(traindataselected {:,:});TrainDataNormalized =(TrainDataselected {::,:} - 含义)./ sdtrain;COEF = PCA(TRAINDATANORMALIZE);据/pre>
PCA1 =(featureSelected {:,:} - meanTrain)./ sdTrain * COEF(:,1);PCA2 =(featureSelected {:,:} - meanTrain)./ sdTrain * COEF(:,2);据/pre>
图numdata = size(featureTable,1);散射(PCA1,PCA2,[],1:NUMDATA,据span style="color:#A020F0">'填充'据/span>)Xlabel(据span style="color:#A020F0">'PCA 1'据/span>)ylabel(据span style="color:#A020F0">'PCA 2'据/span>cbar =彩色杆;ylabel(c bar,[据span style="color:#A020F0">'时间 ('据/span>时髦据span style="color:#A020F0">')'据/span>])据/pre>
HealthIndicator = PCA1;据/pre>
图绘图(PeazerfeLed.date,HealthIndicator,据span style="color:#A020F0">'-O'据/span>)Xlabel(据span style="color:#A020F0">'时间'据/span>) 标题(据span style="color:#A020F0">'健康指标'据/span>)据/pre>
适合剩余使用寿命(RUL)估计的符合指数劣化模型据/h3>
HealthIndicator = HealthIndicator - HealthIndicator(1);据/pre>
阈值=健康指标(结束);据/pre>
MDL = exponentialDegradationModel(据span style="color:#0000FF">......据/span>'theta'据/span>,1,据span style="color:#0000FF">......据/span>'θthetavariance'据/span>,1E6,据span style="color:#0000FF">......据/span>'beta'据/span>,1,据span style="color:#0000FF">......据/span>'etavariance'据/span>,1E6,据span style="color:#0000FF">......据/span>“Phi”据/span>,-1,据span style="color:#0000FF">......据/span>'noisevariance'据/span>,(0.1 *阈值/(阈值+ 1))^ 2,据span style="color:#0000FF">......据/span>'slopedetectionlevel'据/span>,0.05);据/pre>
在每次迭代%保存记录据/span>totalDay =长度(healthIndicator) - 1;estRULs =零(totalDay,1);trueRULs =零(totalDay,1);CIRULs =零(totalDay,2);pdfRULs =细胞(totalDay,1);据span style="color:#228B22">%为绘图更新创建数字和轴据/span>图AX1 =子图(2,1,1);Ax2 =子图(2,1,2);据span style="color:#0000FF">为了据/span>CURRENTDAY = 1:totalDay据span style="color:#228B22">%更新模型参数后验分布据/span>更新(MDL,[CURRENTDAY healthIndicator(CURRENTDAY)])据span style="color:#228B22">%预测剩下的使用寿命据/span>[estrul,cirul,pdfrul] = predictrul(mdl,据span style="color:#0000FF">......据/span>[截止日期医疗器(截止日期)],据span style="color:#0000FF">......据/span>临界点);Truerul = Totaldday - Currentday + 1;据span style="color:#228B22">%更新rul分布图据/span>Helperplottrend(AX1,截止日,Healthindicator,MDL,阈值,时级);Helperplotrul(Ax2,Truerul,Estrul,Cirul,Pdfrul,TimeUnit)据span style="color:#228B22">%保持预测结果据/span>estRULs(currentDay)=estRUL;特鲁尔(当前日)=特鲁尔;CIRULs(当前日期:)=CIRUL;pdfRULs{currentDay}=pdfRUL;据span style="color:#228B22">%暂停0.1秒让动画可见据/span>暂停(0.1)据span style="color:#0000FF">结尾据/span>
绩效分析据/h3>
alpha = 0.2;Detecttime = mdl.slopedetectionstant;prob = helperalphalambdaplot(alpha,trueruls,estruls,ciruls,据span style="color:#0000FF">......据/span>pdfruls,检测时间,断点,时髦);标题(据span style="color:#A020F0">“\α-\拉姆达块”据/span>)据/pre>
图t = 1:overgay;抓住据span style="color:#A020F0">在据/span>绘图(t,prob)图([断点断点],[0 1],据span style="color:#A020F0">'k-。'据/span>) 抓住据span style="color:#A020F0">离开据/span>Xlabel([据span style="color:#A020F0">'时间 ('据/span>时髦据span style="color:#A020F0">')'据/span>])ylabel(据span style="color:#A020F0">'可能性'据/span>) 传奇(据span style="color:#A020F0">“预测rul内的概率”据/span>那据span style="color:#A020F0">“列车试验断点”据/span>) 标题([据span style="color:#A020F0">“\ Alpha绑定中的概率,\ alpha ='据/span>num2str(alpha * 100)据span style="color:#A020F0">'%'据/span>])据/pre>
工具书类据/h3>
exponentialDegradationModel据/code>