風力タービン高速ベアリングの経過予測

この例では、風力タービン ベアリングの残存耐用期間 (RUL) をリアルタイムで予測する指数劣化モデルの作成方法を説明します。指数劣化モデルはそのパラメーターの事前確率分布と最新の測定値に基づいて RUL を予測します (故障に至るまで実行された履歴データはモデル パラメーターの事前確率分布を推定するのに役立ちますが、必須ではありません)。モデルは、重要な劣化トレンドをリアルタイムで検出でき、新しい観測データが利用可能になると、そのパラメーターの事前確率分布を更新します。この例は、一般的な予測のワークフローに従っています。これにはデータのインポートと調査、特徴の抽出と後処理、特徴量重要度のランク付けと融合、モデルの当てはめと予測、および性能の解析が含まれます。

データセット

データセットは 20 歯のピニオン ギアで駆動する 2MW 風力タービンの高速シャフトから収集されたものです [1]。1 日あたり 6 秒間の振動信号を連続 50 日間取得しています (3 月 17 日には 2 つの測定値があり、この例ではこれを 2 日分として扱います)。50 日の期間中に内輪の不具合が発生しベアリングの故障の原因となっています。

ツールボックスではデータセットのコンパクトなバージョンを利用できます。コンパクトなデータセットを使用するには、データセットを現在のフォルダーにコピーしてその書き込み権限を有効にします。

copyfile(...
    fullfile(matlabroot, 'toolbox', 'predmaint', ...
    'predmaintdemos', 'windTurbineHighSpeedBearingPrognosis'), ...
    'WindTurbineHighSpeedBearingPrognosis-Data-main')
fileattrib(fullfile('WindTurbineHighSpeedBearingPrognosis-Data-main', '*.mat'), '+w')

コンパクトなデータセットの測定タイム ステップは 5 日です。

timeUnit = '\times 5 day';

完全なデータセットを使う場合は、https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data のリンクからリポジトリ全体を zip ファイルとしてダウンロードし、このライブ スクリプトと同じディレクトリに保存します。次のコマンドを使用してファイルを解凍します。完全なデータセットの測定タイム ステップは 1 日です。

if exist('WindTurbineHighSpeedBearingPrognosis-Data-main.zip', 'file')
    unzip('WindTurbineHighSpeedBearingPrognosis-Data-main.zip')
    timeUnit = 'day';
end

この例の結果は完全なデータセットから生成されています。この例を実行するには完全なデータセットをダウンロードすることを強く推奨します。コンパクトなデータセットでは有意義な結果を得られない可能性があります。

データのインポート

風力タービン データの fileEnsembleDatastore を作成します。データには振動信号とタコメーター信号が含まれています。fileEnsembleDatastore がファイル名を解析して日付情報を IndependentVariables として抽出します。詳細については、この例に関連する補助ファイルにある補助関数を参照してください。

hsbearing = fileEnsembleDatastore(...
    fullfile('.', 'WindTurbineHighSpeedBearingPrognosis-Data-main'), ...
    '.mat');
hsbearing.DataVariables = ["vibration", "tach"];
hsbearing.IndependentVariables = "Date";
hsbearing.SelectedVariables = ["Date", "vibration", "tach"];
hsbearing.ReadFcn = @helperReadData;
hsbearing.WriteToMemberFcn = @helperWriteToHSBearing;
tall(hsbearing)

Starting parallel pool (parpool) using the 'Processes' profile ...
Connected to parallel pool with 4 workers.

ans =

  M×3 tall table

            Date                vibration             tach     
    ____________________    _________________    ______________

    11-Mar-2013 03:00:24    {292968×1 double}    {187×1 double}
    16-Mar-2013 06:56:43    {292968×1 double}    {180×1 double}
    21-Mar-2013 00:33:14    {292968×1 double}    {188×1 double}
    26-Mar-2013 01:41:50    {292968×1 double}    {180×1 double}
    31-Mar-2013 19:38:18    {292968×1 double}    {180×1 double}
    05-Apr-2013 21:35:37    {292968×1 double}    {182×1 double}
    10-Apr-2013 23:14:07    {292968×1 double}    {187×1 double}
    15-Apr-2013 22:55:24    {292968×1 double}    {167×1 double}
             :                      :                  :
             :                      :                  :

振動信号のサンプル レートは 97656 Hz です。

fs = 97656; % Hz

データの調査

この節では、データを時間領域と周波数領域の両方で調査して、予測をするためにどの特徴を抽出するかのヒントを探します。

まず、振動信号を時間領域で可視化します。このデータセットには、50 日間連続で測定された 6 秒間の振動信号が 50 個あります。そこで、50 の振動信号を順次プロットします。

reset(hsbearing)
tstart = 0;
figure
hold on
while hasdata(hsbearing)
    data = read(hsbearing);
    v = data.vibration{1};
    t = tstart + (1:length(v))/fs;
    % Downsample the signal to reduce memory usage
    plot(t(1:10:end), v(1:10:end))
    tstart = t(end);
end
hold off
xlabel('Time (s), 6 second per day, 50 days in total')
ylabel('Acceleration (g)')

時間領域の振動信号から、信号のインパルス性が増加のトレンドにあることがわかります。尖度、ピークツーピーク値、クレスト ファクターなどの信号のインパルス性を定量化するインジケーターは、この風力タービン ベアリング データセットの今後の見通しに役立つ特徴の候補となります [2]。

一方でスペクトル尖度は、周波数領域で風力タービンの今後の見通しを立てる強力なツールと見なされます [3]。スペクトル尖度の変化を時間に沿って可視化するには、スペクトル尖度の値を周波数と測定日の関数としてプロットします。

hsbearing.DataVariables = ["vibration", "tach", "SpectralKurtosis"];
colors = parula(50);
figure
hold on
reset(hsbearing)
day = 1;
while hasdata(hsbearing)
    data = read(hsbearing);
    data2add = table;
    
    % Get vibration signal and measurement date
    v = data.vibration{1};
    
    % Compute spectral kurtosis with window size = 128
    wc = 128;
    [SK, F] = pkurtosis(v, fs, wc);
    data2add.SpectralKurtosis = {table(F, SK)};
    
    % Plot the spectral kurtosis
    plot3(F, day*ones(size(F)), SK, 'Color', colors(day, :))
    
    % Write spectral kurtosis values
    writeToLastMemberRead(hsbearing, data2add);
    
    % Increment the number of days
    day = day + 1;
end
hold off
xlabel('Frequency (Hz)')
ylabel('Time (day)')
zlabel('Spectral Kurtosis')
grid on
view(-45, 30)
cbar = colorbar;
ylabel(cbar, 'Fault Severity (0 - healthy, 1 - faulty)')

カラー バーに示される故障重大度では、測定の日付が 0 ～ 1 のスケールに正規化されています。マシンの状態が劣化するにつれ、スペクトル尖度の値は 10 kHz 付近で徐々に増加することが観察されます。平均や標準偏差などのスペクトル尖度の統計的特徴は、ベアリング劣化のインジケーターの候補となります [3]。

特徴抽出

前節での解析に基づいて、時間領域の信号およびスペクトル尖度から求めた統計的特徴のコレクションを抽出します。特徴に関する数学的な詳細については、[2-3] を参照してください。

まず、特徴の名前を、fileEnsembleDatastore に書き込む前に DataVariables にあらかじめ割り当てます。

hsbearing.DataVariables = [hsbearing.DataVariables; ...
    "Mean"; "Std"; "Skewness"; "Kurtosis"; "Peak2Peak"; ...
    "RMS"; "CrestFactor"; "ShapeFactor"; "ImpulseFactor"; "MarginFactor"; "Energy"; ...
    "SKMean"; "SKStd"; "SKSkewness"; "SKKurtosis"];

各アンサンブル メンバーについて特徴の値を計算します。

hsbearing.SelectedVariables = ["vibration", "SpectralKurtosis"];
reset(hsbearing)
while hasdata(hsbearing)
    data = read(hsbearing);
    v = data.vibration{1};
    SK = data.SpectralKurtosis{1}.SK;
    features = table;
    
    % Time Domain Features
    features.Mean = mean(v);
    features.Std = std(v);
    features.Skewness = skewness(v);
    features.Kurtosis = kurtosis(v);
    features.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);
    
    % Spectral Kurtosis related features
    features.SKMean = mean(SK);
    features.SKStd = std(SK);
    features.SKSkewness = skewness(SK);
    features.SKKurtosis = kurtosis(SK);
    
    % write the derived features to the corresponding file
    writeToLastMemberRead(hsbearing, features);
end

独立変数 Date と、すべての抽出された特徴を選択して、特徴テーブルを作成します。

hsbearing.SelectedVariables = ["Date", "Mean", "Std", "Skewness", "Kurtosis", "Peak2Peak", ...
    "RMS", "CrestFactor", "ShapeFactor", "ImpulseFactor", "MarginFactor", "Energy", ...
    "SKMean", "SKStd", "SKSkewness", "SKKurtosis"];

特徴テーブルは小さく (50 行 15 列) メモリに十分収まるので、処理の前にテーブルを集めます。ビッグ データの場合は、出力がメモリに十分収まることが確実になるまで、tall 形式で操作を実行することをお勧めします。

featureTable = gather(tall(hsbearing));

Evaluating tall expression using the Parallel Pool 'Processes':
- Pass 1 of 1: Completed in 5.2 sec
Evaluation completed in 5.2 sec

table を timetable に変換して、時間情報が常に特徴の値に関連付けられるようにします。

featureTable = table2timetable(featureTable)

featureTable=10×15 timetable
            Date             Mean       Std       Skewness     Kurtosis    Peak2Peak     RMS      CrestFactor    ShapeFactor    ImpulseFactor    MarginFactor      Energy        SKMean       SKStd      SKSkewness    SKKurtosis
    ____________________    _______    ______    __________    ________    _________    ______    ___________    ___________    _____________    ____________    __________    __________    ________    __________    __________

    11-Mar-2013 03:00:24     0.2139     2.089     0.0065791     3.0405      21.217      2.0999      4.9387         1.2556          6.2009           3.7076       1.2919e+06      0.011681    0.042011      -0.7614        7.566  
    16-Mar-2013 06:56:43    0.23281    1.9755    -0.0060687     3.0069      17.336      1.9892      4.3183          1.254          5.4151           3.4137       1.1592e+06     0.0081655    0.040512       1.0274       6.4863  
    21-Mar-2013 00:33:14    0.18899    2.1852      0.000367     3.1416      24.884      2.1934      6.0332         1.2592          7.5972           4.3616       1.4094e+06    0.00085468    0.066465     -0.38397       11.257  
    26-Mar-2013 01:41:50    0.25741    2.2293     0.0042305     3.0975      23.712      2.2441      5.3639         1.2575          6.7451           3.7797       1.4754e+06      0.011485    0.040831      0.17366       3.3943  
    31-Mar-2013 19:38:18    0.25027    2.1337     -0.003749     3.0971      20.513      2.1483      5.1751         1.2583          6.5119           3.8141       1.3521e+06      0.015608    0.045412       1.5794       7.4012  
    05-Apr-2013 21:35:37    0.21185    2.2492     0.0060094     3.3807      25.088      2.2592      5.7628         1.2691          7.3138           4.1087       1.4953e+06      0.047117     0.12901       3.0512       13.269  
    10-Apr-2013 23:14:07    0.33425    2.6118     0.0021607     3.8872      33.833      2.6331      6.6096         1.2837          8.4847           4.1365       2.0312e+06      0.061578     0.19793       3.7348       17.826  
    15-Apr-2013 22:55:24    0.35205    2.0334     -0.011765      3.938      26.445      2.0636      6.4372         1.2869          8.2839           5.1658       1.2476e+06      0.068291     0.20724       2.9265       10.459  
    20-Apr-2013 15:13:07    0.15898    2.4311     -0.010162     4.6055      33.622      2.4363      7.2259          1.303          9.4153           5.0356       1.7389e+06       0.11731     0.35557       3.0585       11.716  
    25-Apr-2013 23:22:02    0.25785    2.9787      0.025931      5.437      43.445      2.9899      7.6824         1.3298          10.216           4.5439       2.6189e+06       0.16564     0.52757       3.5712       16.232  

特徴の後処理

抽出した特徴には、通常はノイズが関連付けられています。逆のトレンドをもつノイズは、RUL の予測に悪影響を与えることがあります。また、次に導入される、特徴の性能メトリクスの 1 つである単調性は、ノイズに対しロバストでありません。したがって、5 ステップの遅延ウィンドウをもつ因果移動平均フィルターを、抽出した特徴に適用します。ここで「因果」とは、移動平均フィルター処理に将来の値を使用しないという意味です。

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

以下は平滑化の前と後の特徴を示す例です。

figure
hold on
plot(featureTable.Date, featureTable.SKMean)
plot(featureTableSmooth.Date, featureTableSmooth.SKMean)
hold off
xlabel('Time')
ylabel('Feature Value')
legend('Before smoothing', 'After smoothing')
title('SKMean')

移動平均の平滑化によって信号に時間遅延が生じますが、RUL の予測で適切なしきい値を選択することによって遅延の効果を緩和できます。

学習データ

実際には、今後の見通しのアルゴリズムを開発している最中にライフサイクル全体のデータを利用することはできません。しかし、ライフサイクルの初期段階でいくらかのデータが収集済みであると仮定するのは妥当です。したがって、最初の 20 日間 (ライフサイクルの 40%) に収集されたデータを学習データとして扱います。次に示す特徴量重要度のランク付けと融合は、学習データのみに基づいています。

breaktime = datetime(2013, 3, 27);
breakpoint = find(featureTableSmooth.Date < breaktime, 1, 'last');
trainData = featureTableSmooth(1:breakpoint, :);

特徴量重要度のランク付け

この例では、[3] に提案されている単調性を使用して、今後の見通しに用いる特徴の有用性を定量化します。

$\mathit{i}$ 番目の特徴 ${\mathit{x}}_{\mathit{i}}$ の単調性は次のように計算します。

$$\mathrm{Monotonicity}\left({\mathit{x}}_{\mathit{i}}\right)=\frac{1}{\mathit{m}}\sum _{\mathit{j}=1}^{\mathit{m}}\frac{|\mathrm{number}\mathrm{of}\mathrm{positive}\mathrm{diff}\left({\mathit{x}}_{\mathit{i}}^{\mathit{j}}\right)-\mathrm{number}\mathrm{of}\mathrm{negative}\mathrm{diff}\left({\mathit{x}}_{\mathit{i}}^{\mathit{j}}\right)|}{\mathit{n}-1}$$

$\mathit{n}$ は測定点の数で、ここでは $\mathit{n}=50$ です。$\mathit{m}$ は監視されるマシンの数で、ここでは $\mathit{m}=1$ です。${\mathit{x}}_{\mathit{i}}^{\mathit{j}}$ は $\mathit{j}$ 番目のマシンで測定された $\mathit{i}$ 番目の特徴です。$\mathrm{diff}\left({\mathit{x}}_{\mathit{i}}^{\mathit{j}}\right)={\mathit{x}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}\right)-{\mathit{x}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}-1\right)$ であり、これはすなわち信号 ${\mathit{x}}_{\mathit{i}}^{\mathit{j}}$ の差分です。

% Since moving window smoothing is already done, set 'WindowSize' to 0 to
% turn off the smoothing within the function
featureImportance = monotonicity(trainData, 'WindowSize', 0);
helperSortedBarPlot(featureImportance, 'Monotonicity');

信号の尖度が、単調性に基づいた最上位の特徴です。

特徴量重要度スコアが 0.3 より大きい特徴は、次節での特徴融合の対象として選択されます。

trainDataSelected = trainData(:, featureImportance{:,:}>0.3);
featureSelected = featureTableSmooth(:, featureImportance{:,:}>0.3)

featureSelected=10×15 timetable
            Date             Mean       Std       Skewness     Kurtosis    Peak2Peak     RMS      CrestFactor    ShapeFactor    ImpulseFactor    MarginFactor      Energy       SKMean       SKStd      SKSkewness    SKKurtosis
    ____________________    _______    ______    __________    ________    _________    ______    ___________    ___________    _____________    ____________    __________    _________    ________    __________    __________

    11-Mar-2013 03:00:24     0.2139     2.089     0.0065791     3.0405      21.217      2.0999      4.9387         1.2556          6.2009           3.7076       1.2919e+06     0.011681    0.042011      -0.7614        7.566  
    16-Mar-2013 06:56:43    0.22336    2.0322    0.00025522     3.0237      19.277      2.0445      4.6285         1.2548           5.808           3.5607       1.2255e+06    0.0099232    0.041261        0.133       7.0262  
    21-Mar-2013 00:33:14     0.2119    2.0832    0.00029248      3.063      21.146      2.0941      5.0967         1.2563          6.4044           3.8276       1.2868e+06    0.0069003    0.049662    -0.039322       8.4363  
    26-Mar-2013 01:41:50    0.22328    2.1197      0.001277     3.0716      21.787      2.1316      5.1635         1.2566          6.4896           3.8157        1.334e+06    0.0080464    0.047455     0.013923       7.1758  
    31-Mar-2013 19:38:18    0.22868    2.1225    0.00027179     3.0767      21.532       2.135      5.1658         1.2569           6.494           3.8153       1.3376e+06    0.0095588    0.047046      0.32701       7.2209  
    05-Apr-2013 21:35:37    0.22587    2.1437     0.0012281     3.1274      22.125      2.1557      5.2653          1.259          6.6306           3.8642       1.3639e+06     0.015819    0.060706      0.78105       8.2289  
    10-Apr-2013 23:14:07    0.24593    2.2308    0.00049166     3.2685      24.228      2.2445      5.5438         1.2636          7.0113           3.9357       1.4871e+06     0.024135    0.086693       1.5304       9.9389  
    15-Apr-2013 22:55:24     0.2658    2.2404    -0.0004577     3.4237      25.746       2.257      5.8969         1.2691          7.4894           4.2277       1.5018e+06     0.034156     0.11448       1.8469       10.601  
    20-Apr-2013 15:13:07     0.2608    2.2814    -0.0022125     3.6677      27.202      2.2974      6.0957         1.2764          7.7924           4.3401       1.5568e+06     0.053565     0.16267       2.4207       10.677  
    25-Apr-2013 23:22:02    0.26087    2.4063     0.0014042     4.0576      30.491      2.4217      6.4821         1.2885          8.3709           4.4674       1.7473e+06     0.079258     0.24379       2.9869       12.817  

次元削減と特徴の融合

この例では、次元削減と特徴の融合に主成分分析 (PCA) を使用します。PCA を実行する前に、特徴を同じスケールに正規化することをお勧めします。PCA の係数と、正規化に使用される平均および標準偏差は学習データから取得され、データセット全体に適用されます。

meanTrain = mean(trainDataSelected{:,:});
sdTrain = std(trainDataSelected{:,:});
trainDataNormalized = (trainDataSelected{:,:} - meanTrain)./sdTrain;
coef = pca(trainDataNormalized);

平均、標準偏差、および PCA 係数は、データセット全体を処理するために使用されます。

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

データを最初の 2 つの主成分の空間で可視化します。

figure
numData = size(featureTable, 1);
scatter(PCA1, PCA2, [], 1:numData, 'filled')
xlabel('PCA 1')
ylabel('PCA 2')
cbar = colorbar;
ylabel(cbar, ['Time (' timeUnit ')'])

プロットは、マシンが故障に近づくにつれて第 1 主成分が増加することを示しています。したがって、第 1 主成分は融合された健康インジケーターの有望な候補となります。

healthIndicator = PCA1;

健康インジケーターを可視化します。

figure
plot(featureSelected.Date, healthIndicator, '-o')
xlabel('Time')
title('Health Indicator')

残存耐用期間 (RUL) 推定の指数劣化モデルの当てはめ

指数劣化モデルは次のように定義されます。

$$\mathit{h}\left(\mathit{t}\right)=\varphi +\theta \exp \left(\beta \mathit{t}+ϵ-\frac{{\sigma }^2}{2}\right)$$

ここで、$\mathit{h}\left(\mathit{t}\right)$ は時間の関数として表した健康インジケーターです。$\varphi$ は、定数と見なされる切片項です。$\theta$ と $\beta$ はモデルの傾きを決定する乱数パラメーターで、$\theta$ は対数正規分布、$\beta$ はガウス分布です。各タイム ステップ $\mathit{t}$ で、$\theta$ および $\beta$ の分布は $\mathit{h}\left(\mathit{t}\right)$ の最新の観測値に基づいて事後に更新されます。$ϵ$ はガウス ホワイト ノイズで、$\mathit{N}\left(0,{\sigma }^2\right)$ に従います。指数関数の $-\frac{{\sigma }^2}{2}$ 項は、$\mathit{h}\left(\mathit{t}\right)$ の期待値が次を満たすようにするものです。

$\mathit{E}\left[\mathit{h}\left(\mathit{t}\right)|\theta ,\beta \right]=\varphi +\theta \exp \left(\beta \mathit{t}\right)$.

ここで指数劣化モデルを前節で抽出した健康インジケーターに当てはめ、次の節で性能を評価します。

まず、健康インジケーターをシフトさせて 0 から開始するようにします。

healthIndicator = healthIndicator - healthIndicator(1);

しきい値の選択は、通常はマシンの過去の記録または何らかのドメイン固有の知識に基づいています。このデータセットでは履歴データが利用できないので、健康インジケーターの最後の値をしきい値として選択します。平滑化の遅延効果が一部緩和されるよう、平滑化された (履歴) データに基づいたしきい値を選択することをお勧めします。

threshold = healthIndicator(end);

履歴データを利用できる場合は、exponentialDegradationModel が提供する fit メソッドを使用して事前確率分布と切片を推定します。しかし、この風力タービン ベアリングのデータセットでは履歴データが利用できません。傾きパラメーターの事前確率分布は大きな分散で任意に選択されるので ($\mathit{E}\left(\theta \right)=1,\mathrm{Var}\left(\theta \right)={10}^6,\mathit{E}\left(\beta \right)=1,\mathrm{Var}\left(\beta \right)={10}^6$)、モデルは主に観測されたデータに依存しています。$\mathit{E}\left[\mathit{h}\left(0\right)\right]=\varphi +\mathit{E}\left(\theta \right)$ に基づき、切片 $\varphi$ を $-1$ に設定して、モデルも 0 から開始されるようにします。

健康インジケーターの変化とノイズの変化間の関係は、次のように求めることができます。

$$\Delta \mathit{h}\left(\mathit{t}\right)\approx \left(\mathit{h}\left(\mathit{t}\right)-\varphi \right)\Delta ϵ\left(\mathit{t}\right)$$

ここで、しきい値の付近では、ノイズの標準偏差が健康インジケーターの変化の 10% の原因となると仮定します。したがって、ノイズの標準偏差を $\frac{10\%\cdot \mathrm{threshold}}{\mathrm{threshold}-\varphi }$ と表すことができます。

指数劣化モデルは傾きの重要性を評価する機能も提供します。健康インジケーターの重要な傾きが検出されると、モデルは前の観測値を破棄し、元の事前確率分布に基づいて推定を再び開始します。検出アルゴリズムの感度は、SlopeDetectionLevel を指定することで調整できます。p 値が SlopeDetectionLevel 未満の場合、傾きが検出されたとして宣言されます。ここで SlopeDetectionLevel は 0.05 に設定されます。

次に、上で述べたパラメーターを使って指数劣化モデルを作成します。

mdl = exponentialDegradationModel(...
    'Theta', 1, ...
    'ThetaVariance', 1e6, ...
    'Beta', 1, ...
    'BetaVariance', 1e6, ...
    'Phi', -1, ...
    'NoiseVariance', (0.1*threshold/(threshold + 1))^2, ...
    'SlopeDetectionLevel', 0.05);

predictRUL メソッドと update メソッドを使用して RUL を予測し、パラメーター分布をリアルタイムで更新します。

% Keep records at each iteration
totalDay = length(healthIndicator) - 1;
estRULs = zeros(totalDay, 1);
trueRULs = zeros(totalDay, 1);
CIRULs = zeros(totalDay, 2);
pdfRULs = cell(totalDay, 1);

% Create figures and axes for plot updating
figure
ax1 = subplot(2, 1, 1);
ax2 = subplot(2, 1, 2);

for currentDay = 1:totalDay
    
    % Update model parameter posterior distribution
    update(mdl, [currentDay healthIndicator(currentDay)])
    
    % Predict Remaining Useful Life
    [estRUL, CIRUL, pdfRUL] = predictRUL(mdl, ...
                                         [currentDay healthIndicator(currentDay)], ...
                                         threshold);
    trueRUL = totalDay - currentDay + 1;
    
    % Updating RUL distribution plot
    helperPlotTrend(ax1, currentDay, healthIndicator, mdl, threshold, timeUnit);
    helperPlotRUL(ax2, trueRUL, estRUL, CIRUL, pdfRUL, timeUnit)
    
    % Keep prediction results
    estRULs(currentDay) = estRUL;
    trueRULs(currentDay) = trueRUL;
    CIRULs(currentDay, :) = CIRUL;
    pdfRULs{currentDay} = pdfRUL;
    
    % Pause 0.1 seconds to make the animation visible
    pause(0.1)
end

Warning: Data does not satisfy exponential model assumption. Check model's prior parameters or try a different degradation model.

Warning: Data does not satisfy exponential model assumption. Check model's prior parameters or try a different degradation model.

以下はリアルタイムの RUL 推定のアニメーションです。

性能解析

$\alpha$-$\lambda$ プロットが今後の見通しの性能解析に使用されます [5]。ここで、$\alpha$ の範囲は 20% に設定されます。推定される RUL が真の RUL の $\alpha$ の範囲内にある確率は、モデルの性能メトリクスとして次のように計算されます。

$$\mathrm{Pr}\left({\mathit{r}}^{*}\left(\mathit{t}\right)-\alpha {\mathit{r}}^{*}\left(\mathit{t}\right)<\mathit{r}\left(\mathit{t}\right)<{\mathit{r}}^{*}\left(\mathit{t}\right)+\alpha {\mathit{r}}^{*}\left(\mathit{t}\right)|\Theta \left(\mathit{t}\right)\right)$$

ここで、$\mathit{r}\left(\mathit{t}\right)$ は時間 $\mathit{t}$ での推定された RUL、${\mathit{r}}^{*}\left(\mathit{t}\right)$ は時間 $\mathit{t}$ での真の RUL、$\Theta \left(\mathit{t}\right)$ は時間 $\mathit{t}$ での推定されたモデル パラメーターです。

alpha = 0.2;
detectTime = mdl.SlopeDetectionInstant;
prob = helperAlphaLambdaPlot(alpha, trueRULs, estRULs, CIRULs, ...
    pdfRULs, detectTime, breakpoint, timeUnit);
title('\alpha-\lambda Plot')

事前設定された事前確率分布は真の事前確率分布を反映しないため、モデルは通常、適切なパラメーター分布に調整されるまで数タイム ステップを要します。利用できるデータ点が増えるほど予測は正確になります。

予測された RUL が $\alpha$ の範囲内にある確率を可視化します。

figure
t = 1:totalDay;
hold on
plot(t, prob)
plot([breakpoint breakpoint], [0 1], 'k-.')
hold off
xlabel(['Time (' timeUnit ')'])
ylabel('Probability')
legend('Probability of predicted RUL within \alpha bound', 'Train-Test Breakpoint')
title(['Probability within \alpha bound, \alpha = ' num2str(alpha*100) '%'])

参考文献

[1] Bechhoefer, Eric, Brandon Van Hecke, and David He."Processing for improved spectral analysis."Annual Conference of the Prognostics and Health Management Society, New Orleans, LA, Oct. 2013.

[2] Ali, Jaouher Ben, et al."Online automatic diagnosis of wind turbine bearings progressive degradations under real experimental conditions based on unsupervised machine learning."Applied Acoustics 132 (2018):167-181.

[3] Saidi, Lotfi, et al."Wind turbine high-speed shaft bearings health prognosis through a spectral Kurtosis-derived indices and SVR."Applied Acoustics 120 (2017):1-8.

[4] Coble, Jamie Baalis."Merging data sources to predict remaining useful life–an automated method to identify prognostic parameters." (2010).

[5] Saxena, Abhinav, et al."Metrics for offline evaluation of prognostic performance."International Journal of Prognostics and Health Management 1.1 (2010):4-23.

参考

exponentialDegradationModel

トピック

• 残存耐用期間の予測のための特徴選択
• RUL 推定器モデルを使用した RUL の推定
• 類似性ベースの残存耐用期間推定