风力涡轮机高速承载预后

打开直播脚本

该实例说明了如何建立指数退化模型来实时预测风力机轴承的剩余使用寿命。指数退化模型基于其参数先验和最新的度量来预测RUL(历史运行到故障数据可以帮助估计模型参数先验，但它们不是必需的)。该模型能够实时检测显著的退化趋势，并在有新的观测数据时更新其参数先验值。该示例采用了典型的预测流程:数据导入和探索、特征提取和后处理、特征重要性排序和融合、模型拟合和预测，以及性能分析。

数据集

数据集是从2MW风力涡轮机高速轴由20齿小齿轮[1]驱动收集。6秒的振动信号被获取的每一天为50个连续天（有2个测量3月17日，其在该示例中视为两个天）。内圈故障开发和跨越导致50天期间，轴承的故障。

Toolbox中提供了一个紧凑的数据集版。要使用Compact DataSet，将DataSet复制到当前文件夹并启用其写入权限。

copyfile（..．fullfile (matlabroot'工具箱'，“predmaint”，..．'predmaintdemos'，“windTurbineHighSpeedBearingPrognosis”），..．'WindTurbineHighSpeedBearingPrognosis-DATA-硕士fileattrib（fullfile（'WindTurbineHighSpeedBearingPrognosis-DATA-硕士，'* .mat'），' + w '）

Compact DataSet的测量时间步长为5天。

TimeUnit =.'\ times 5天'；

对于整个数据集，请访问以下链接https://github.com/mathworks/windturbinehighspeedbearingprognosiss-data.，将整个存储库作为zip文件下载，并将其保存在与此实时脚本相同的目录中。使用此命令解压缩文件。完整数据集的测量时间步长为1天。

如果存在（'windturbinehighpeedbearingprognosis-data-master.zip'，“文件”）解压缩（'windturbinehighpeedbearingprognosis-data-master.zip') timeUnit ='天'；结束

本例中的结果是从完整的数据集生成的。强烈建议下载完整的数据集来运行此示例。从紧凑数据集生成的结果可能没有意义。

数据导入

创建一个fileEnsembleDatastore风力发电机的数据。数据包括一个振动信号和一个转速表信号。的fileEnsembleDatastore将解析文件名并将日期信息提取为IndependentVariables．有关更多细节，请参阅与本示例相关的支持文件中的helper函数。金宝app

hsbearing = fileEnsembleDatastore (..．fullfile (“。”，'WindTurbineHighSpeedBearingPrognosis-DATA-硕士），..．'.mat'）;hsbearing.datavariables = [“振动”，“tach”];hsbearing。IndependentVariables =“日期”；hsbearing.SelectedVariables = [“日期”，“振动”，“tach”];hsbearing.readfcn = @helperreaddata;hsbearing.writetomemberfcn = @helperwriteTohsbearing;高大（Hsbearing）

ans = M×3高表日期振动性心动过速  ____________________ _________________ _______________ 07 - 3月- 2013 01:57:46(585936×1双)(2446×1双)08 - mar - 2013 02:34:21(585936×1双)(2411×1双)09 - 3月- 2013 02:33:43(585936×1双)(2465×1双)10 - 3月- 2013年03:01:02(585936×1双)(2461×1双)11 - 3月- 2013 03:00:24[585936×1[2506×1 double] 12-Mar-2013 06:17:10 [585936×1 double] [2447×1 double] 13-Mar-2013 06:34:04 [585936×1 double] [2438×1 double] 14-Mar-2013 06:50:41 [585936×1 double] [2390×1 double]::::::

振动信号的采样率为97656 Hz。

FS = 97656;％Hz.

数据探索

本节探讨时域和频域的数据，并寻找用于预测目的的特征提取的灵感。

首先在时域中可视化振动信号。在该数据集中，在连续50天中有50个振动信号6秒。现在绘制50振动信号彼此一个。

重置（hsbearing）tstart = 0;图持有在而Hasdata (hsbearing) data = read(hsbearing);v = data.vibration {1};t = tstart +（1：长度（v））/ fs;向下采样信号以减少内存使用Plot (t(1:10:end)， v(1:10:end)) tstart = t(end);结束持有关闭Xlabel（“时间(秒)，每天6秒，共50天”）ylabel（“加速度(g)”）

时域振动信号表现出信号冲动性增大的趋势。量化信号冲动性的指标，如峰度，峰峰值，波峰因子等等，对于该风力涡轮机轴承数据集进行潜在的预后特征[2]。

另一方面，谱峰度被认为是在频域[3]进行风力机预测的有力工具。为了直观地看到谱峰度随时间的变化，绘制谱峰度值作为频率和测量日的函数。

hsbearing.datavariables = [“振动”，“tach”，“SpectralKurtosis”];颜色= parula（50）;图持有在复位（hsbearing）天= 1;而Hasdata (hsbearing) data = read(hsbearing);data2add =表;％获得振动信号和测量日期v = data.vibration {1};%计算谱峰度，窗宽= 128WC = 128;[SK，F] = pkurtosis（V，FS，WC）;data2add.SpectralKurtosis = {表（F，SK）};％绘制光谱峰度Plot3（F，Day * Ones（尺寸（f）），sk，'颜色'，颜色（日，:)）%写谱峰度值WriteTolastmemberRead（Hsbearing，Data2Add）;%增加天数Day = Day + 1;结束持有关闭Xlabel（'频率（Hz）'）ylabel（的时间(天)）Zlabel（“谱峰态”）网格在查看（-45,30）cbar =彩色杆;ylabel（c bar，“故障级别(0 -正常，1 -故障)”）

颜色条中显示的故障严重程度是归一化为0到1刻度的测量日期。观察到，随着机器条件的退化，光谱峰度值在10 kHz左右逐渐增加。谱峰度的统计特征，如平均值、标准差等等，将是轴承劣化的潜在指标[3]。

特征提取

基于前一节中的分析，将提取源自时域信号和光谱峰度的统计特征的集合。[2-3]中提供了有关该功能的更多数学细节。

首先，在将其写入FileSeneMbleyataStore之前，在DataVaria中预先分配特征名称。

hsbearing.datavariables = [hsbearing.datavariables;..．“意思”；“std”；“偏态”；“峰度”；“peak2peak”；..．“rms”；“crestfactor”；“shapefactor”；“ImpulseFactor”；“marginfactor”；“能量”；..．“skmean”；“skstd”；“skskewness”；“skkurtosis”];

计算每个集合成员的特征值。

hsbearing.SelectedVariables = [“振动”，“SpectralKurtosis”];重置(hsbearing)而Hasdata (hsbearing) data = read(hsbearing);v = data.vibration {1};SK = data.SpectralKurtosis {1} .SK;特点=表;%时域特性特征.. 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);光谱峰度相关特征特征.skmean =平均（sk）;特征.skstd = std（sk）;特征.skskewness = skewness（sk）;特征.skkurtosis = kurtosis（sk）;%将派生的特性写入相应的文件writeToLastMemberRead (hsbearing、特点);结束

选择独立变量日期并将提取的所有特征构造特征表。

hsbearing.SelectedVariables = [“日期”，“意思”，“std”，“偏态”，“峰度”，“peak2peak”，..．“rms”，“crestfactor”，“shapefactor”，“ImpulseFactor”，“marginfactor”，“能量”，..．“skmean”，“skstd”，“skskewness”，“skkurtosis”];

由于要素表足够小以适合内存（50到15），因此在处理前收集表。对于大数据，建议以高表格执行操作，直到您相信输出足够小以适合内存。

featureTable =收集(高(hsbearing));

使用并行池“本地”评估高表达： -  PASS 1为1：1秒评估完成1秒

将表转换为时间表，以便时间信息始终与特征值相关联。

featureTable = table2timetable（featureTable）

featureTable =50×15时间表日均值标准偏度峰度Peak2Peak RMS峰值因子形状因子ImpulseFactor MarginFactor能源SKMean SKStd SKSkewness SKKurtosis ____________________ _______ ______ ___________ ________ _________ ______ ___________ ___________ _____________ ____________ __________ __________ ________ __________ __________ 07月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二点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-三月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时50分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点50分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 6点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 18时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 0点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分五十秒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⋮

特征后处理

提取的特征通常与噪声相关联。具有相反趋势的噪音有时可能对rul预测有害。另外，要引入的特征性能指标，单调性待引入的一个是不稳健的噪声。因此，将具有5个步骤的滞后窗口的因果移动平均滤波器应用于提取的特征，其中“因果”表示移动平均滤波中没有将未来的值用于未来的值。

variableNames = featureTable.Properties.VariableNames;featureTableSmooth = varfun(@(x) movmean(x， [5 0])， featureTable);featureTableSmooth.Properties.VariableNames = variableNames;

这是一个示例，示出了平滑之前和之后的功能。

图持有在图（featureTable.Date，featureTable.SKMean）图（featureTableSmooth.Date，featureTableSmooth.SKMean）保持关闭Xlabel（“时间”）ylabel（“特征值”)传说(在平滑的，'平滑后'）标题（“SKMean”）

移动平均平滑引入了信号的时间延迟，但是通过在RUL预测中选择适当的阈值，可以减轻延迟效果。

培训数据

在实践中，在开发预测算法时，无法获得整个生命周期的数据，但可以合理地假设已经收集了生命周期早期的一些数据。因此，前20天(生命周期的40%)收集的数据被视为训练数据。下面的特征重要性排序和融合仅基于训练数据。

breaktime = datetime（2013,3,27）;breakpoint = find（featureTablesmooth.date “最后一次”）;TrainData = FeatureTablesSmooth（1：断点，:);

功能重要性排名

在这个例子中，[3]提出的单调性被用来量化特征对预后的价值。

单调性的 $我$ 特征 $x_{我}$ 是计算

$单调性（ x_{我} ）＝ \frac{1}{米} \sum_{j ＝ 1}^{米} \frac{| 号码 of 积极 diff （ x_{我}^{j} ） - 号码 of 否定 diff （ x_{我}^{j} ） |}{n - 1}$

在哪里 $n$ 是测量点数，在这种情况下 $n ＝ 50$ ． $米$ 在这种情况下，监控的机器数量是多少 $米＝ 1$ ． $x_{我}^{j}$ 是 $我$ 测量的功能 $j$ 机。 $diff （ x_{我}^{j} ）＝ x_{我}^{j} （ t ） - x_{我}^{j} （ t - 1 ）$ ，即信号的差异 $x_{我}^{j}$ ．

由于移动窗口平滑已经完成，因此将“Windowsize”设置为0到0%关闭函数内的平滑featureImportance =单调性(trainData,'Windowsize'，0）;helperSortedBarPlot（featureImportance，的单调性）;

信号的峰度是信号单调性的最主要特征。

选择具有大于0.3的特征重要性分数的功能对于下一节中的功能融合。

TrainDataselected = TrainData（:, FeedingImportance {：，：}> 0.3）;feationElted = featureTablesmooth（:, FeedingImportance {：，：}> 0.3）

feedingeled =50×5时间表日均值峭度形状因子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六点五十分41秒0.23299 3.0249 1.2544 3.9047 0.02824 15-MAR-2013六时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 0时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月2013年0点39分五十○秒0.23079 3.1481 1.2593 4.072 0.058908⋮

降维与特征融合

本例中采用主成分分析(PCA)进行降维和特征融合。在进行PCA之前，将特征归一化到相同的尺度是一个很好的做法。注意，归一化使用的PCA系数和均值、标准差都是从训练数据中得到的，并应用于整个数据集。

含义=卑鄙（TrainDataselected {：，：，：}）;sdtrain = std（traindataselected {：，：}）;TrainDataNormalized =（TrainDataselected {：：，：}  - 含义）./ sdtrain;COEF = PCA（TRAINDATANORMALIZE）;

平均值，标准偏差和PCA系数用于处理整个数据集。

PCA1 = (featureSelected{:，:} - meanTrain) ./ sdTrain * coef(:， 1);PCA2 = (featureSelected{:，:} - meanTrain) ./ sdTrain * coef(:， 2);

可视化前两个主组件的空间中的数据。

图numdata = size（featureTable，1）;散射（PCA1，PCA2，[]，1：NUMDATA，“填充”)包含('PCA 1'）ylabel（《PCA 2》cbar =彩色杆;ylabel（c bar，['时间（'timeUnit'）']）

图中显示，当机器接近故障时，第一个主成分是增加的。因此，第一个主成分是一个很有前途的融合健康指标。

HealthIndicator = PCA1;

可视化健康指标。

图绘图（PeazerfeLed.date，HealthIndicator，'-O')包含(“时间”）标题（'健康指标'）

适合剩余使用寿命（RUL）估计的符合指数劣化模型

指数劣化模型被定义为

$h （ t ）＝ ϕ + θ. 经验值（ β t + ϵ - \frac{{σ.}^{2}}{2} ）$

在哪里 $h （ t ）$ 是健康指标作为时间的函数。 $ϕ$ 被认为是一个常数截距项。 $θ.$ 和 $β$ 是否随机参数决定了模型的斜率，在哪里 $θ.$ 是Lognormal分布的 $β$ 是系统。在每个时间步 $t$ ，分布 $θ.$ 和 $β$ 是根据最近的观察更新到后部的 $h （ t ）$ ． $ϵ$ 高斯白噪声是否屈服于 $N （ 0 ， {σ.}^{2} ）$ ．的 $- \frac{{σ.}^{2}}{2}$ 指数中的术语是预期 $h （ t ）$ 满足

$E ［ h （ t ） | θ. ， β ］＝ ϕ + θ. 经验值（ β t ）$ ．

这里，指数劣化模型适合最后一节中提取的健康指示符，并且在下一节中评估性能。

首先移动运行状况指示器，使其从0开始。

HealthIndicator = HealthIndicator  -  HealthIndicator（1）;

阈值的选择通常基于机器的历史记录或某种域特定知识。由于此数据集中未使用历史数据，因此选择健康指示符的最后值作为阈值。建议基于平滑（历史）数据选择阈值，以便部分减轻平滑的延迟效果。

阈值= HealthIndicator（END）;

如果有历史数据，请使用适合提供的方法exponentialDegradationModel估计先验和截距。然而，历史数据不能用于该风力涡轮机轴承数据集。斜率参数的先验是任意选取的，方差很大( $E （ θ. ）＝ 1 ， var. （ θ. ）＝ 10^{6} ， E （ β ）＝ 1 ， var. （ β ）＝ 10^{6}$ ）因此该模型主要依赖于观察到的数据。基于 $E ［ h （ 0 ）］＝ ϕ + E （ θ. ）$ ，拦截 $ϕ$ 设置为 $- 1$ 这样模型也会从0开始。

健康指标变异与噪声变异之间的关系可以推导为

$δ. h （ t ） \approx （ h （ t ） - ϕ ） δ. ϵ （ t ）$

这里假设噪声的标准偏差在接近阈值时导致健康指示器的10％变化。因此，噪声的标准偏差可以表示为 $\frac{10 ％ \cdot 门槛}{门槛 - ϕ}$ ．

指数退化模型还提供了评估斜坡重要性的功能。一旦检测到健康指标的显著斜率，模型将忘记以前的观察结果，并基于原始先验重新开始估计。检测算法的灵敏度可以通过指定来调整SlopeDetectionLevel．如果P值小于SlopeDetectionLevel，则宣布探测到坡度。在这里SlopeDetectionLevel设置为0.05。

现在使用上面讨论的参数创建指数劣化模型。

mdl = exponentialDegradationModel (..．“θ”, 1..．'θthetavariance'，1e6，..．“β”, 1..．'etavariance'，1e6，..．“φ”，-1，..．'noisevariance'，（0.1 *阈值/（阈值+ 1））^ 2，..．'slopedetectionlevel'，0.05）;

使用predictRUL和更新方法预测RUL并实时更新参数分布。

在每次迭代％保存记录totalDay = length(healthIndicator) - 1;estRULs = 0 (totalDay, 1);trueRULs = 0 (totalDay, 1);CIRULs = zeros(totalDay, 2);pdfRULs = cell(totalDay, 1); / /累计天数%创建图形和坐标轴用于绘图更新图AX1 =子图（2,1,1）;Ax2 =子图（2,1,2）;对于currentDay = 1: totalDay%更新模型参数后验分布更新（MDL，[CURRENTDAY healthIndicator（CURRENTDAY）]）%预测剩余的使用寿命[estRUL, CIRUL, pdfRUL] = prestrul (mdl, pdl)..．[截止日期医疗器（截止日期）]，..．阈值);true = totalDay - currentDay + 1;％更新rul分布图helperplotrend (ax1, currentDay, healthIndicator, mdl, threshold, timeUnit);helperPlotRUL(ax2, trueRUL, estRUL, CIRUL, pdfRUL, timeUnit)保持预测结果雌激流（截止日期）= estrul;Trueruls（截止日期）= Truerul;线圈（截止日，:) = cirul;pdfruls {curstday} = pdfrul;％暂停0.1秒让动画可见暂停(0.1)结束

这里是实时RUL估计的动画。

性能分析

$α.$ - $λ.$ 情节用于预后性能分析[5]，其中 $α.$ Bound设置为20%。估计的RUL在 $α.$ RUL的边界作为模型的性能度量来计算:

$公关（ r^{＊} （ t ） - α. r^{＊} （ t ） < r （ t ） < r^{＊} （ t ） + α. r^{＊} （ t ） | θ. （ t ））$

在哪里 $r （ t ）$ 是估计的统治 $t$ ， $r^{＊} （ t ）$ 是真正的rul $t$ ， $θ. （ t ）$ 是估计的模型参数吗？ $t$ ．

α= 0.2;detectTime = mdl.SlopeDetectionInstant;prob = helperAlphaLambdaPlot(alpha, trueRULs, estRULs, CIRULs，..．pdfruls，检测时间，断点，时髦）;标题（“α- \ \λ阴谋”）

由于预设事先不反映真实的先前，模型通常需要几个时间步骤来调整到适当的参数分布。随着更多数据点可用，预测变得更加准确。

将预测的RUL的概率可视化 $α.$ 绑定。

图t = 1：overgay;持有在绘图（t，prob）图（[断点断点]，[0 1]，'k-。'抱紧关闭包含(['时间（'timeUnit'）']）ylabel（'概率')传说(“在\alpha范围内预测RUL的概率”，“Train-Test断点”）标题（['概率在\ α范围内，\ α = 'num2str（alpha * 100）'％']）