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

この例では、風力タービン ベアリングの残存耐用期間 (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)

ans =

  M×3 tall table

            Date                vibration             tach      
    ____________________    _________________    _______________

    07-Mar-2013 01:57:46    {585936×1 double}    {2446×1 double}
    08-Mar-2013 02:34:21    {585936×1 double}    {2411×1 double}
    09-Mar-2013 02:33:43    {585936×1 double}    {2465×1 double}
    10-Mar-2013 03:01:02    {585936×1 double}    {2461×1 double}
    11-Mar-2013 03:00:24    {585936×1 double}    {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 日間連続で測定された 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 2 sec
Evaluation completed in 2 sec

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

featureTable = table2timetable(featureTable)

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

    07-Mar-2013 01:57:46    0.34605    2.2705      0.0038699     2.9956      21.621      2.2967      4.9147         1.2535          6.1607           3.3625       3.0907e+06     0.0013007     0.02575     -0.22011       3.3392  
    08-Mar-2013 02: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.0046258    0.020869     0.057756       3.3494  
    09-Mar-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.6208e+06    -0.0010934    0.022712     -0.49972       4.9732  
    10-Mar-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.3894e+06     0.0087136    0.034486       1.4755       8.1605  
    11-Mar-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.515e+06       0.01358    0.032145     0.096699       3.8647  
    12-Mar-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.5833e+06     0.0017442    0.028323     -0.13925       3.8028  
    13-Mar-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.3051e+06     0.0038714    0.029292       -1.476       8.1441  
    14-Mar-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.9575e+06     0.0013091     0.02249      0.27383       2.8652  
    15-Mar-2013 06: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.0023288    0.047793      -2.6038       20.273  
    16-Mar-2013 06: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.0075703     0.03041      0.51657       3.9987  
    17-Mar-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.6437e+06     0.0084583    0.037896       2.3773       11.521  
    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.019554    0.055275       3.1725       17.624  
    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.016271    0.064586       2.8774       11.658  
    20-Mar-2013 00: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.0011064    0.051602     -0.06065       7.0731  
    21-Mar-2013 00: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.0039722     0.06632     -0.40713       12.154  
    22-Mar-2013 00:39:50    0.20569    2.1861      0.0020375     3.2691      26.439      2.1958      6.1753         1.2633          7.8011           4.4882        2.825e+06      0.020589    0.077706       2.5994       11.084  
      ⋮

特徴の後処理

抽出した特徴には、通常はノイズが関連付けられています。逆のトレンドをもつノイズは、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=50×5 timetable
            Date             Mean      Kurtosis    ShapeFactor    MarginFactor     SKStd  
    ____________________    _______    ________    ___________    ____________    ________

    07-Mar-2013 01:57:46    0.34605     2.9956       1.2535          3.3625        0.02575
    08-Mar-2013 02:34:21    0.29507     3.0075        1.254          3.5428       0.023309
    09-Mar-2013 02:33:43    0.26962     3.0125        1.254          3.6541        0.02311
    10-Mar-2013 03:01:02    0.25565     3.0197       1.2544          3.7722       0.025954
    11-Mar-2013 03:00:24    0.24756     3.0247       1.2546          3.7793       0.027192
    12-Mar-2013 06:17:10    0.25519     3.0236       1.2544          3.7479       0.027381
    13-Mar-2013 06:34:04      0.233     3.0272       1.2545          3.8282       0.027971
    14-Mar-2013 06:50:41    0.23299     3.0249       1.2544          3.9047       0.028241
    15-Mar-2013 06:50:03     0.2315      3.033       1.2548          3.9706       0.032421
    16-Mar-2013 06:56:43    0.23475     3.0273       1.2546          3.9451       0.031742
    17-Mar-2013 06:56:04    0.23498     3.0407       1.2551          3.9924         0.0327
    17-Mar-2013 18:47:56    0.21839     3.0533       1.2557          3.9792       0.037192
    18-Mar-2013 18:47:15    0.21943     3.0778       1.2567          4.0538       0.043075
    20-Mar-2013 00:33:54    0.23854     3.0905       1.2573          3.9658       0.047927
    21-Mar-2013 00:33:14    0.23537     3.1044       1.2578          3.9862       0.051015
    22-Mar-2013 00:39:50    0.23079     3.1481       1.2593           4.072       0.058897
      ⋮

次元削減と特徴の融合

この例では、次元削減と特徴の融合に主成分分析 (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

以下はリアルタイムの 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 の推定
• 類似性ベースの残存耐用期間推定