风力发电机高速轴承预测- MATLAB和Simulink万博1manbetx - 万博1manbetx,s manbetx 845,万博尤文图斯

风力发电机高速轴承预后

开立真实脚本

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

数据集

数据集收集自一个2MW的风力涡轮机高速轴驱动的20齿小齿轮[1]。连续50天每天采集6秒的振动信号(3月17日有2次测量，在本例中为2天)。在50天的时间内，轴承发生了内圈故障，导致轴承失效。

工具箱中提供了数据集的压缩版本。若要使用压缩数据集，请将数据集复制到当前文件夹并启用其写权限。

拷贝文件（...fullfile (matlabroot'工具箱',“predmaint”,...“predmaintdemos”,“windTurbineHighSpeedBearingPrognosis”）...“WindTurbineHighSpeedBearingPrognosis-Data-master”)fileattrib (fullfile (“WindTurbineHighSpeedBearingPrognosis-Data-master”,'*。垫'）' + w ')

紧凑数据集的测量时间步长为5天。

timeUnit =“\ * 5天”;

对于完整的数据集，请转到这个链接https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data，下载整个仓库作为一个zip文件，并在同一目录下保存为这个活的脚本。解压缩使用此命令的文件。对全部数据集测定时间的步骤是1天。

如果存在('WindTurbineHighSpeedBearingPrognosis-Data-master.zip',“文件”)解压缩('WindTurbineHighSpeedBearingPrognosis-Data-master.zip')timeUnit =“天”;结束

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

数据导入

创建一个fileEnsembleDatastore风力涡轮机的数据。数据包含一个振动信号和一个转速计信号。的fileEnsembleDatastore将解析的文件名并提取日期信息IndependentVariables.有关详细信息，请参阅与此示例关联的支持文件中的帮助函数。万博1manbetx

hsbearing = fileEnsembleDatastore (...fullfile (“。”,“WindTurbineHighSpeedBearingPrognosis-Data-master”）...'。垫');hsbearing.DataVariables = [“振动”,“一环”];hsbearing。IndependentVariables =“日期”;hsbearing。选择edVariables = [“日期”,“振动”,“一环”];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双)12 - 3月- 2013年06:17:10(585936×1双)(2447×1双)13 - 3月- 2013 06:34:04(585936×1双)(2438×1双)14 - 3月- 201306:50:41[585936×1双][2390×1双]:::::::

振动信号采样率为97656 Hz。

fs = 97656;%赫兹

数据探索

本节将探索时域和频域的数据，并从中寻找对预后有帮助的特征。

第一可视化在时域的振动的信号。在该数据集中，存在50个连续天测量6秒50个的振动信号。现在绘制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);结束持有离包含(时间，每天6秒，总共50天）ylabel（“加速度(g)”)

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

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

hsbearing.DataVariables = [“振动”,“一环”,“SpectralKurtosis”];颜色= parula（50）;身材保持上复位(hsbearing)天= 1;而hasdata(hsbearing) data = read(hsbearing);data2add =表;％获取振动信号和测量日期V = data.vibration {1};计算窗口大小为128的光谱峰度wc = 128;[SK, F] =峰度(v, fs, wc);data2add。{table(F, SK)};%表示光谱峰度plot3 (F,天*的(大小(F)), SK,“颜色”颜色(天,:))写谱峰度值writeToLastMemberRead（hsbearing，data2add）;%增加天数day = day + 1;结束持有离包含(的频率(赫兹)）ylabel（的时间(天)）zlabel（“谱峰态”网格)上视图(- 45,30)cbar = colorbar;ylabel (cbar“故障严重程度(0 -健康，1 -故障)”)

颜色栏中显示的故障严重程度是归一化为0到1刻度的测量数据。可以看出，随着机器状态的下降，10khz左右的谱峰度值逐渐增大。光谱峰度的统计特征，如均值、标准差等等等，将是轴承退化的潜在指标[3]。

特征提取

在前一节分析的基础上，提取时域信号和光谱峰度的统计特征集合。关于特征的更多数学细节在[2-3]中提供。

首先，在将其写入fileEnsembleDatastore之前，在DataVariables中预先分配特性名称。

hsbearing.DataVariables = [hsbearing.DataVariables;...“的意思是”;“性病”;“偏态”;“峰度”;“Peak2Peak”;...“RMS”;“峰值因数”;“ShapeFactor”;“ImpulseFactor”;“MarginFactor”;“能量”;...“SKMean”;“SKStd”;“SKSkewness”;“SKKurtosis”];

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

hsbearing。选择edVariables = [“振动”,“SpectralKurtosis”];重置(hsbearing)而hasdata(hsbearing) data = read(hsbearing);v = data.vibration {1};SK = data.SpectralKurtosis {1} .SK;特点=表;%时域特征特性。的意思= (v);特性。Std =性病(v);特性。偏态=偏斜度(v);特性。峰度=峰度(v);特性。Peak2Peak = Peak2Peak (v); features.RMS = rms(v); features.CrestFactor = max(v)/features.RMS; features.ShapeFactor = features.RMS/mean(abs(v)); features.ImpulseFactor = max(v)/mean(abs(v)); features.MarginFactor = max(v)/mean(abs(v))^2; features.Energy = sum(v.^2);%光谱峰度相关特征features.SKMean =平均值（SK）;features.SKStd = STD（SK）;features.SKSkewness =偏斜度（SK）;features.SKKurtosis =峭度（SK）;将派生的特性写入对应的文件writeToLastMemberRead (hsbearing、特点);结束

选择自变量日期和所有提取的特征，构造特征表。

hsbearing。选择edVariables = [“日期”,“的意思是”,“性病”,“偏态”,“峰度”,“Peak2Peak”,...“RMS”,“峰值因数”,“ShapeFactor”,“ImpulseFactor”,“MarginFactor”,“能量”,...“SKMean”,“SKStd”,“SKSkewness”,“SKKurtosis”];

由于特征表是小到足以装配在存储器（50 15），则处理前收集表。对于大数据，建议在高大的格式，直到你有信心，输出足够小，以适应在内存中进行操作。

featureTable =收集(高(hsbearing));

评估使用并行池“本地”高表达： - 的1遍1：在1秒完成评价在1秒完成

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

featureTable = table2timetable（featureTable）

featureTable =50×15的时间表日期的意思是性病偏态峰度Peak2Peak RMS CrestFactor ShapeFactor ImpulseFactor MarginFactor能源SKMean SKStd SKSkewness SKKurtosis ____________________ _________…………………………_________________ ________ ____ _______ __________ __________ ________ __________ __________ 07 - 3月- 2013 01:57:46 0.34605 2.2705 0.0038699 2.9956 21.621 2.2967 0.001253 4.9147 1.2535 6.1607 3.3625 3.0907 e + 06 - mar - 2013 02:34:21 08年0.025674 -0.22882 3.362 0.24409 2.06210.0030103 3.0195 19.31 2.0765 4.9129 1.2545 6.163 3.7231 2.5266 e + 06年0.0046823 0.020888 0.057651 3.3508 09 - 3月- 2013 02:33:43 0.21873 2.1036 -0.0010289 3.0224 21.474 2.1149 5.2143 1.2539 6.5384 3.8766 2.6208 e + 06 -0.0011084 0.022705 -0.50004 4.9953 10 - 3月- 2013年03:01:02 0.21372 2.0081 0.001477 3.0415 19.52 2.0194 5.286 1.2556 6.637 4.1266 2.3894 e + 06 0.0087035 0.034456 1.4705 8.1235 11 - 3月- 2013 03:00:24 0.21518 2.0606 0.0010116 3.0445 21.217 2.0718 5.0058 1.2554 6.2841 3.8078 2.515 0.013559 e + 060.032193 0.11355 3.848 12 - 3月- 2013年06:17:10 0.29335 2.0791 -0.008428 3.018 20.05 2.0997 4.7966 1.2537 6.0136 3.5907 2.5833 e + 06年0.0017539 0.028326 -0.1288 3.8072 13 - 3月- 2013 06:34:04 0.21293 1.972 -0.0014294 3.0174 18.837 1.9834 4.8496 1.2539 6.0808 3.8441 2.3051 e + 06 0.0039353 - 0.029292 -1.4734 - 8.1242 14 - 3月- 2013 06:50:41 0.24401 1.8114 0.0022161 3.0057 17.862 1.8278 4.8638 1.2536 6.0975 4.1821 1.9575 e + 06 15 - 3月- 2013 06:50:03 0.0013107 0.022468 0.27438 2.8597 0.20984 1.9973 0.001559 3.0711 - 21.120.0023769 2.0083 5.4323 1.2568 6.8272 4.2724 2.3633 e + 06 16 - 3月- 2013 06:56:43 0.047767 -2.5832 20.171 0.23318 1.9842 -0.0019594 3.0072 18.832 1.9979 5.0483 1.254 6.3304 3.9734 2.3387 e + 06年0.0076327 0.030418 0.52322 4.0082 17 - 3月- 2013 06:56:04 0.21657 2.113 -0.0013711 3.1247 21.858 2.1241 5.4857 1.2587 6.9048 4.0916 2.6437 e + 06年0.0084907 0.037876 2.3753 11.491 17 - 3月- 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.6891 e + 06年0.019624 0.055537 3.1986 17.79618 - 3月- 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.6824 e + 20 - 3月- 2013年06 0.016315 0.064516 2.8735 11.632 00:33:54 0.35865 2.2536 -0.002308 3.0817 22.633 2.2819 0.0011097 5.2771 1.2571 6.6339 3.6546 3.0511 e + 06 21 - 3月- 2013 00:33:14 0.051539 -0.056774 7.0712 0.1908 2.1782 -0.00019286 3.1548 25.515 2.1865 6.0521 1.26 7.6258 4.3945 2.8013 0.0040392 e + 06 22 - 3月- 2013 00:39:50 0.066254 -0.39587 12.111 0.20569 2.1861 0.0020375 3.2691 26.439 2.1958 6.17531.2633 7.8011 4.4882 2.825 e + 06年0.020676 0.077728 2.6038 11.088⋮

特征后处理

提取的特征通常与噪声相关。趋势相反的噪声有时会对预测结果产生不利影响。另外，接下来要介绍的特征性能指标之一单调性对噪声不具有鲁棒性。因此，我们对提取的特征使用了一个延迟窗口为5步的因果移动均值滤波，其中“因果”表示在移动均值滤波中没有使用未来值。

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

下面是一个显示平滑前后特征的例子。

身材保持上情节(featureTable。日期,featureTable.SKMean) plot(featureTableSmooth.Date, featureTableSmooth.SKMean) hold离包含(“时间”）ylabel（“特征值”)传说(在平滑的,平滑后的)标题('SKMean')

移动平均的平滑介绍了信号的时间延迟，但该延迟的效果可通过在RUL预测选择合适的阈值而减轻。

训练数据

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

时间=日期时间(2013年3月27日);断点=找到(featureTableSmooth。日期<休息时间，1，“最后一次”);trainData = featureTableSmooth（1：断点，:);

功能重要性排名

在本例中，我们使用[3]提出的单调性来量化这些特征的优点，以达到预测的目的。

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

$单调性 (x_{我}) = \frac{1}{米} •_{j = 1}^{米} \frac{| 数的正 diff (x_{我}^{j}) - 数的负 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到关闭函数内的平滑featureImportance =单调性(trainData,“WindowSize”,0);helperSortedBarPlot (featureImportance的单调性);

信号的峰度是基于单调性的最上面的特征。

与特征重要性得分特征大于0.3被选择用于在下一节特征融合。

trainDataSelected = trainData（:, featureImportance {：，：}> 0.3）;featureSelected = featureTableSmooth（:, featureImportance {：，：}> 0.3）

featureSelected =50×5时间表日期的意思是峰度ShapeFactor MarginFactor SKStd ____________________ _________ _______ ________ ________ ___________ 07 - 3月- 2013 01:57:46 08年0.34605 2.9956 1.2535 3.3625 0.025674 - 3月- 2013 02:34:21 0.29507 3.0075 1.254 3.5428 0.023281 09 - 3月- 2013 02:33:43 0.26962 3.0125 1.254 3.6541 0.023089 10 - 3月- 2013年11 03:01:02 0.25565 3.0197 1.2544 3.7722 0.025931 - 3月- 2013 03:00:24 0.24756 3.0247 1.2546 3.7793 0.027183 3月12 - 13 - 2013 06:17:10 0.25519 3.0236 1.2544 3.7479 0.027374 0.233 - 3月- 2013 06:34:043.0272 1.2545 3.8282 0.027977 14-Mar-2013 06:50:41 0.23299 3.0249 1.2544 3.9047 0.02824 15-Mar-2013 06:50:03 0.2315 3.033 1.2548 3.9706 0.032417 16-Mar-2013 06:56:43 0.23475 3.0273 1.2546 3.9451 0.031744 17-Mar-2013 06: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 00:33:54 0.23854 3.0905 1.2573 3.9658 0.047942 21-Mar-2013 00:33:14 0.23537 3.1044 1.2578 3.9862 0.051023 22-Mar-2013 00:39:50 0.23079 3.1481 1.2593 4.072 0.058908 ⋮

降维与特征融合

在本例中，我们使用主成分分析(PCA)进行降维和特征融合。在执行PCA之前，最好将特征归一化到相同的范围内。注意，PCA系数和归一化中使用的均值和标准差是从训练数据中获得的，并应用于整个数据集。

meanTrain =平均值（trainDataSelected {：，：}）;sdTrain = STD（trainDataSelected {：，：}）;trainDataNormalized =（trainDataSelected {：，：}  -  meanTrain）./ sdTrain;COEF = PCA（trainDataNormalized）;

使用均值、标准差和主成分分析系数来处理整个数据集。

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

在前两个主成分的空间可视化的数据。

图numData = SIZE（featureTable，1）;散射（PCA1，PCA2，[]，1：numData，“填充”）xlabel（“PCA 1”）ylabel（《PCA 2》) cbar = colorbar;ylabel (cbar ['时间 （'TIMEUNIT“)”])

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

healthIndicator = PCA1;

可视化的健康指标。

图绘制(featureSelected。日期、healthIndicator'-o'）xlabel（“时间”)标题(“健康指示器”)

对于剩余使用寿命（RUL）估计拟合指数降解模型

指数退化模型定义为

$h (t) = ϕ + θ 经验值 (β t + ϵ - \frac{σ^{2}}{2})$

在哪里 $h (t)$ 是健康指标，是时间的函数。 $ϕ$ 是被认为是常数的截距项。 $θ$ 和 $β$ 模型的斜率是由随机参数决定的吗 $θ$ lognormal-distributed和 $β$ 是系统。在每个时间步长 $t$ ，的分布 $θ$ 和 $β$ 是根据最新的观测结果更新到后面的吗 $h (t)$ . $ϵ$ 高斯白噪声是否屈服于 $N (0, σ^{2})$ .的 $- \frac{σ^{2}}{2}$ 指数项的期望值是 $h (t)$ 满足

$E (h (t) | θ, β] = ϕ + θ 经验值 (β t)$ .

这里，指数退化模型适合于上一节中提取的健康指标，并在下一节中评估其性能。

首先移动健康指示器，使它从0开始。

healthIndicator = healthIndicator - healthIndicator(1);

阈值的选择通常是基于机器或某些特定领域知识的历史记录。由于没有历史数据在该数据集可用，健康指示符的最后值被选择为阈值。它建议选择基于平滑（历史的）数据，以使得平滑的延迟效果将被部分减轻所述阈值。

阈值= healthIndicator(结束);

如果历史数据可用，就使用它适合方法提供了一种通过exponentialDegradationModel估计先验和截距。然而，历史数据不能用于这个风力发电机轴承数据集。斜率参数的先验值是随机选取的，其方差较大( $E (θ) = 1, Var (θ) = 10^{6}, E (β) = 1, Var (β) = 10^{6}$ ），使得模型大多依赖于所观察到的数据。基于 $E (h (0)] = ϕ + E (θ)$ ,拦截 $ϕ$ 被设置为 $- 1$ 所以模型也会从0开始。

健康指示符的变化和噪声的变化之间的关系可被推导为

$Δ h (t) • (h (t) - ϕ) Δ ϵ (t)$

这里假定的噪声的标准偏差，以使健康指示符的变化的10％时，它是接近阈值。因此，噪声的标准偏差可以被表示为 $\frac{10 % \cdot 阈值}{阈值 - ϕ}$ .

指数退化模型还提供了一种评估斜坡重要性的功能。一旦检测到健康指标的显著斜率，模型就会忘记之前的观测结果，重新基于原始先验进行估计。检测算法的灵敏度可以通过指定进行调整SlopeDetectionLevel.如果p值小于SlopeDetectionLevel，则表示斜坡已被探测。在这里SlopeDetectionLevel设置为0.05。

现在创建上面讨论的参数的指数衰减模型。

mdl = exponentialDegradationModel (...“θ”,1...“ThetaVariance”1 e6,...“β”,1...“BetaVariance”1 e6,...“φ”,1...“NoiseVariance”，（0.1 *阈值/（阈值+ 1））^ 2，...“SlopeDetectionLevel”,0.05);

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

在每次迭代中保持记录totalDay =长度(healthIndicator) - 1;estRULs = 0 (totalDay, 1);trueRULs = 0 (totalDay, 1);CIRULs = 0 (totalDay, 2);pdfRULs = cell(totalDay, 1);为情节的更新创建图形和轴图AX1 =副区（2,1，1）;AX2 =副区（2,1，2）;为currentDay = 1: totalDay更新模型参数后验分布更新(mdl [currentDay healthIndicator (currentDay)))%预测剩余使用寿命[estRUL, CIRUL, pdfRUL] = predictRUL(mdl，...[CURRENTDAY healthIndicator（CURRENTDAY）]...阈值);trueRUL = totalDay - currentDay + 1;更新的RUL分布图helperPlotTrend(ax1, currentDay, healthIndicator, mdl, threshold, timeUnit);helperPlotRUL(ax2, trueRUL, estRUL, CIRUL, pdfRUL, timeUnit)%保持预测结果estRULs (currentDay) = estRUL;trueRULs (currentDay) = trueRUL;CIRULs(currentDay，:) = CIRUL;pdfRULs {currentDay} = pdfRUL;%暂停0.1秒使动画可见暂停(0.1)结束

下面是实时RUL估计的动画。

性能分析

$α$ - $λ$ plot用于预后性能分析[5]，其中 $α$ 绑定设置为20%。估计的RUL在。之间的概率 $α$ 真实RUL的界计算为模型的性能指标:

$公关 (r^{*} (t) - α r^{*} (t) < r (t) < r^{*} (t) + α r^{*} (t) | Θ (t))$

在哪里 $r (t)$ 估计的RUL是在什么时候 $t$ , $r^{*} (t)$ 真正的规则是什么 $t$ , $Θ (t)$ 是在时间的估计的模型参数 $t .$

α= 0.2;detectTime = mdl.SlopeDetectionInstant;b = helperAlphaLambdaPlot(α、trueRULs、estRULs、CIRULs、...pdfRULs，detectTime，断点，TIMEUNIT）;标题（“α- \ \λ阴谋”)

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

将预测的RUL的概率形象化 $α$ 绑定。

图t = 1:totalDay;持有上情节（T，概率）情节（[断点断点]，[0 1]，“K-”。）保持离包含(['时间 （'TIMEUNIT“)”])ylabel (“概率”)传说(“在\alpha范围内预测RUL的概率”,“Train-Test断点”）标题（['在\alpha范围内的概率，\alpha = 'num2str(α* 100)“%”])