類似性ベースの残存耐用期間推定

ライブスクリプトを開く

この例では、残存耐用期間 (RUL) 推定の完全なワークフローを作成する方法を説明します。これには前処理、トレンドを示す特徴の選択、センサー融合による健康インジケーターの作成、類似性 RUL 推定器の学習、および予知性能の検証の手順が含まれます。ここでは、PHM2008 challenge データセットの学習データを使用します [1]。

データの準備

データセットは小さいので、劣化データ全体をメモリに読み込むことが可能です。データセットをダウンロードして現在のディレクトリに解凍します。補助関数 helperLoadData を使用して学習テキストファイルを読み込み、timetable の cell 配列に変換します。学習データには、故障に至るまで実行された 260 回のシミュレーションのデータが含まれています。この計測値のグループを "アンサンブル" と呼びます。

url = 'https://ssd.mathworks.com/supportfiles/nnet/data/TurbofanEngineDegradationSimulationData.zip';
websave('TurbofanEngineDegradationSimulationData.zip',url);
unzip('TurbofanEngineDegradationSimulationData.zip')
degradationData = helperLoadData('train_FD002.txt');
degradationData(1:5)

ans=5×1 cell array
    {149×26 table}
    {269×26 table}
    {206×26 table}
    {235×26 table}
    {154×26 table}

各アンサンブルメンバーは 26 列の table です。列にはマシン ID、タイムスタンプ、3 つの動作状態、および 21 のセンサー測定値のデータが含まれます。

head(degradationData{1})

    id    time    op_setting_1    op_setting_2    op_setting_3    sensor_1    sensor_2    sensor_3    sensor_4    sensor_5    sensor_6    sensor_7    sensor_8    sensor_9    sensor_10    sensor_11    sensor_12    sensor_13    sensor_14    sensor_15    sensor_16    sensor_17    sensor_18    sensor_19    sensor_20    sensor_21
    __    ____    ____________    ____________    ____________    ________    ________    ________    ________    ________    ________    ________    ________    ________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________

    1      1         34.998            0.84           100          449.44      555.32      1358.6      1137.2       5.48           8       194.64      2222.7      8341.9       1.02         42.02       183.06       2387.7       8048.6       9.3461        0.02          334         2223           100        14.73       8.8071  
    1      2         41.998          0.8408           100             445       549.9      1353.2      1125.8       3.91        5.71       138.51      2211.6        8304       1.02          42.2       130.42       2387.7       8072.3       9.3774        0.02          330         2212           100        10.41       6.2665  
    1      3         24.999          0.6218            60          462.54      537.31      1256.8      1047.5       7.05        9.02       175.71      1915.1      8001.4       0.94         36.69       164.22         2028       7864.9       10.894        0.02          309         1915         84.93        14.08       8.6723  
    1      4         42.008          0.8416           100             445      549.51        1354      1126.4       3.91        5.71       138.46      2211.6        8304       1.02         41.96       130.72       2387.6       8068.7       9.3528        0.02          329         2212           100        10.59       6.4701  
    1      5             25          0.6203            60          462.54      537.07      1257.7      1047.9       7.05        9.03       175.05      1915.1      7993.2       0.94         36.89       164.31         2028       7861.2       10.896        0.02          309         1915         84.93        14.13       8.5286  
    1      6         25.005          0.6205            60          462.54      537.02      1266.4      1048.7       7.05        9.03       175.17      1915.2      7996.1       0.94         36.78       164.27         2028       7868.9       10.891        0.02          306         1915         84.93        14.28        8.559  
    1      7         42.004          0.8409           100             445      549.74      1347.5      1127.2       3.91        5.71       138.71      2211.6      8307.8       1.02         42.19       130.49       2387.7       8075.5       9.3753        0.02          330         2212           100        10.62       6.4227  
    1      8         20.002          0.7002           100          491.19      607.44      1481.7      1252.4       9.35       13.65       334.41      2323.9      8709.1       1.08         44.27       315.11         2388       8049.3       9.2369        0.02          365         2324           100        24.33       14.799

劣化データを、学習データセットと、後で性能評価を行うための検証データセットに分割します。

rng('default')  % To make sure the results are repeatable
numEnsemble = length(degradationData);
numFold = 5;
cv = cvpartition(numEnsemble, 'KFold', numFold);
trainData = degradationData(training(cv, 1));
validationData = degradationData(test(cv, 1));

対象とする変数のグループを指定します。

varNames = string(degradationData{1}.Properties.VariableNames);
timeVariable = varNames{2};
conditionVariables = varNames(3:5);
dataVariables = varNames(6:26);

アンサンブルデータのサンプルを可視化します。

nsample = 10;
figure
helperPlotEnsemble(trainData, timeVariable, ...
    [conditionVariables(1:2) dataVariables(1:2)], nsample)

作業領域クラスタリング

前の節で示したように、故障に至るまで実行された各測定データにおいて劣化プロセスを示す明確なトレンドはありません。この節と次の節では、動作状態を使用してより明確な劣化トレンドをセンサー信号から抽出します。

各アンサンブルメンバーには、"op_setting_1"、"op_setting_2"、"op_setting_3" の 3 つの動作状態が含まれています。まず、各 cell から table を抽出して単一の table に連結します。

trainDataUnwrap = vertcat(trainData{:});
opConditionUnwrap = trainDataUnwrap(:, cellstr(conditionVariables));

すべての操作点を 3 次元散布図に可視化します。明らかに 6 つの領域が示されており、各領域内の点は非常に近接しています。

figure
helperPlotClusters(opConditionUnwrap)

クラスタリング手法を使って 6 つのクラスターを自動的に検出してみましょう。ここでは k-means アルゴリズムを使用します。k-means は最もよく使われるクラスタリングアルゴリズムの 1 つですが、結果が局所的な最適値となる場合があります。k-means クラスタリングアルゴリズムを異なった初期状態で何度か繰り返し、コストが最も低い結果を選択することをお勧めします。ここではアルゴリズムが 5 回実行され、同じ結果となっています。

opts = statset('Display', 'final');
[clusterIndex, centers] = kmeans(table2array(opConditionUnwrap), 6, ...
    'Distance', 'sqeuclidean', 'Replicates', 5, 'Options', opts);

Replicate 1, 1 iterations, total sum of distances = 0.331378.
Replicate 2, 1 iterations, total sum of distances = 0.331378.
Replicate 3, 1 iterations, total sum of distances = 0.331378.
Replicate 4, 1 iterations, total sum of distances = 0.331378.
Replicate 5, 1 iterations, total sum of distances = 0.331378.
Best total sum of distances = 0.331378

クラスタリングの結果と特定されたクラスター重心を可視化します。

figure
helperPlotClusters(opConditionUnwrap, clusterIndex, centers)

プロットに示されているように、クラスタリングアルゴリズムで 6 つの作業領域が正しく検出されます。

作業領域の正規化

異なる作業領域でグループ化された測定値について正規化を実行しましょう。最初に、前の節で検出した作業領域によってグループ化された各センサー測定値の平均と標準偏差を計算します。

centerstats = struct('Mean', table(), 'SD', table());
for v = dataVariables
    centerstats.Mean.(char(v)) = splitapply(@mean, trainDataUnwrap.(char(v)), clusterIndex);
    centerstats.SD.(char(v))   = splitapply(@std,  trainDataUnwrap.(char(v)), clusterIndex);
end
centerstats.Mean

ans=6×21 table
    sensor_1    sensor_2    sensor_3    sensor_4    sensor_5    sensor_6    sensor_7    sensor_8    sensor_9    sensor_10    sensor_11    sensor_12    sensor_13    sensor_14    sensor_15    sensor_16    sensor_17    sensor_18    sensor_19    sensor_20    sensor_21
    ________    ________    ________    ________    ________    ________    ________    ________    ________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________    _________

     489.05      604.91        1502      1311.2      10.52       15.493      394.33        2319      8785.8         1.26      45.487       371.45       2388.2       8135.5       8.6645          0.03      369.69        2319           100       28.527       17.115  
     462.54      536.86      1262.7      1050.2       7.05       9.0277      175.43      1915.4      8015.6      0.93992      36.798       164.57       2028.3       7878.6       10.913          0.02      307.34        1915         84.93       14.263       8.5587  
        445       549.7      1354.4      1127.7       3.91       5.7158      138.63        2212      8328.1       1.0202      42.146       130.55       2388.1       8089.7       9.3745          0.02      331.06        2212           100       10.588       6.3516  
     491.19      607.56      1485.5        1253       9.35       13.657       334.5        2324        8730       1.0777      44.446       314.88       2388.2       8066.1       9.2325      0.022112      365.37        2324           100       24.455       14.674  
     518.67      642.67      1590.3      1408.7      14.62        21.61      553.37      2388.1      9063.3          1.3      47.538       521.41       2388.1         8142       8.4419          0.03      393.18        2388           100        38.82       23.292  
     449.44       555.8      1366.7      1131.5       5.48       8.0003      194.45        2223      8356.4       1.0204      41.981       183.03       2388.1       8072.5       9.3319          0.02      334.22        2223           100       14.831       8.8983

centerstats.SD

ans=6×21 table
     sensor_1     sensor_2    sensor_3    sensor_4     sensor_5     sensor_6     sensor_7    sensor_8    sensor_9    sensor_10     sensor_11    sensor_12    sensor_13    sensor_14    sensor_15    sensor_16     sensor_17    sensor_18    sensor_19     sensor_20    sensor_21
    __________    ________    ________    ________    __________    _________    ________    ________    ________    __________    _________    _________    _________    _________    _________    __________    _________    _________    __________    _________    _________

    1.0346e-11    0.47582      5.8057      8.5359     3.2687e-13    0.0047002    0.67361     0.095286     18.819     0.00017552      0.2538      0.53621     0.098213      16.507      0.038392     6.6965e-16       1.49          0                 0     0.14441      0.08598 
    4.6615e-12    0.36083       5.285      6.8365      1.137e-13    0.0042209    0.44862      0.27006     14.479     0.00087595     0.20843      0.34493      0.28592      13.427      0.043297     3.4003e-16     1.2991          0        3.5246e-12     0.11218     0.066392 
             0    0.43282      5.6422      7.6564     1.1014e-13    0.0049312    0.43972      0.30478     18.325      0.0015563     0.23965      0.34011       0.3286      16.773      0.037476     1.0652e-15     1.4114          0                 0     0.10832     0.064709 
    1.8362e-11    0.46339      5.7888      7.7768     2.5049e-13    0.0047449     0.6046      0.12828      17.58       0.004193     0.23826      0.49843      0.13156      15.457      0.038886       0.004082     1.4729          0                 0     0.13531     0.079592 
    1.2279e-11     0.4967      6.0349      8.9734      6.573e-14    0.0013408    0.86985      0.07082     19.889     2.7314e-14     0.26508      0.73713     0.070545      17.105      0.037289     6.6619e-16     1.5332          0                 0     0.17871       0.1065 
    1.4098e-11     0.4401      5.6013      7.4461     2.2828e-13    0.0017508    0.47564      0.28634     17.446      0.0019595     0.23078      0.37631      0.30766      15.825      0.038139     3.3656e-16     1.4133          0                 0     0.11347     0.067192

各領域の統計値を使用して学習データを正規化することができます。各アンサンブルメンバーについて、各行の操作点を抽出し、その各クラスター中心までの距離を計算して、最も近いクラスター中心を見つけます。次に、各センサー測定値について、平均を減算し、そのクラスターの標準偏差で除算します。標準偏差が 0 に近い場合は、正規化されたセンサー測定値を 0 に設定します。これは、ほぼ一定のセンサー測定値は残存耐用期間の推定に有用でないためです。関数 regimeNormalization の詳細については、最後の節である「補助関数」を参照してください。

trainDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ...
    trainData, 'UniformOutput', false);

作業領域で正規化されたデータを可視化します。一部のセンサー測定値における劣化のトレンドが正規化の後で明らかになっています。

figure
helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(1:4), nsample)

トレンド可能性の解析

次に、トレンドを最もよく示すセンサー測定値を選択して、予測のための健康インジケーターを作成します。各センサー測定値について、線形劣化モデルを推定して信号の勾配をランク付けします。

numSensors = length(dataVariables);
signalSlope = zeros(numSensors, 1);
warn = warning('off');
for ct = 1:numSensors
    tmp = cellfun(@(tbl) tbl(:, cellstr(dataVariables(ct))), trainDataNormalized, 'UniformOutput', false);
    mdl = linearDegradationModel(); % create model
    fit(mdl, tmp); % train mode
    signalSlope(ct) = mdl.Theta;
end
warning(warn);

信号の勾配を並べ替え、勾配が最も大きい 8 つのセンサーを選択します。

[~, idx] = sort(abs(signalSlope), 'descend');
sensorTrended = sort(idx(1:8))

sensorTrended = 8×1

     2
     3
     4
     7
    11
    12
    15
    17

選択された、トレンドを示すセンサー測定値を可視化します。

figure
helperPlotEnsemble(trainDataNormalized, timeVariable, dataVariables(sensorTrended(3:6)), nsample)

トレンドを最もよく示す信号のいくつかは正のトレンドですが、他は負のトレンドである点に注意してください。

健康インジケーターの作成

この節では、センサー測定値を単一の健康インジケーターに融合させる作業に焦点を当てます。インジケーターは類似性ベースのモデルの学習に使用されます。

故障に至るまで実行されたすべてのデータは、健全な状態から開始されると仮定します。初期の健康状態には 1 の値が割り当てられ、故障時の健康状態には 0 の値が割り当てられます。健康状態は時間とともに 1 から 0 へ線形に劣化すると仮定されます。この線形劣化を利用して、センサー測定値の融合を行います。より高度なセンサー融合手法については、文献 [2-5] で説明されています。

for j=1:numel(trainDataNormalized)
    data = trainDataNormalized{j};
    rul = max(data.time)-data.time;
    data.health_condition = rul / max(rul);
    trainDataNormalized{j} = data;
end

健康状態を可視化します。

figure
helperPlotEnsemble(trainDataNormalized, timeVariable, "health_condition", nsample)

すべてのアンサンブルメンバーの健康状態が、さまざまな劣化速度で 1 から 0 に変化します。

次に、Health Condition の線形回帰モデルを、トレンドを最もよく示すセンサー測定値にリグレッサーとして当てはめます。

Health Condition ~ 1 + Sensor2 + Sensor3 + Sensor4 + Sensor7 + Sensor11 + Sensor12 + Sensor15 + Sensor17

trainDataNormalizedUnwrap = vertcat(trainDataNormalized{:});

sensorToFuse = dataVariables(sensorTrended);
X = trainDataNormalizedUnwrap{:, cellstr(sensorToFuse)};
y = trainDataNormalizedUnwrap.health_condition;
regModel = fitlm(X,y);
bias = regModel.Coefficients.Estimate(1)

bias = 0.5000

weights = regModel.Coefficients.Estimate(2:end)

weights = 8×1

   -0.0323
   -0.0300
   -0.0527
    0.0057
   -0.0646
    0.0054
   -0.0431
   -0.0379

センサー測定値をその関連する重みで乗算し、単一の健康インジケーターを作成します。

trainDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), trainDataNormalized, ...
    'UniformOutput', false);

融合された健康インジケーターを学習データについて可視化します。

figure
helperPlotEnsemble(trainDataFused, [], 1, nsample)
xlabel('Time')
ylabel('Health Indicator')
title('Training Data')

複数のセンサーからのデータが単一の健康インジケーターに融合されます。健康インジケーターは、移動平均フィルターで平滑化されています。詳細については、最後の節にある補助関数 "degradationSensorFusion" を参照してください。

検証データへの同じ操作の適用

領域の正規化とセンサー融合のプロセスを検証データセットで繰り返します。

validationDataNormalized = cellfun(@(data) regimeNormalization(data, centers, centerstats), ...
    validationData, 'UniformOutput', false);
validationDataFused = cellfun(@(data) degradationSensorFusion(data, sensorToFuse, weights), ...
    validationDataNormalized, 'UniformOutput', false);

健康インジケーターを検証データについて可視化します。

figure
helperPlotEnsemble(validationDataFused, [], 1, nsample)
xlabel('Time')
ylabel('Health Indicator')
title('Validation Data')

類似性 RUL モデルの作成

次に、学習データを使用して残差ベースの類似性 RUL モデルを作成します。この設定では、モデルは各融合データを 2 次多項式に当てはめようとします。

データ $i$ とデータ $j$ 間の距離は、残差の 1 ノルムにより計算されます。

$d (i, j) = | | y_{j} - {y_{}^{ˆ}}_{j, i} | |_{1}$

ここで、 $y_{j}$ はマシン $j$ の健康インジケーターであり、 ${y_{}^{ˆ}}_{j, i}$ はマシン $i$ で同定された 2 次多項式モデルを使用する、マシン $j$ の推定された健康インジケーターです。

類似性スコアは次の式で計算されます。

$s c o r e (i, j) = e x p (- d (i, j)^{2})$

検証データセット内の 1 つのアンサンブルメンバーについて、モデルは学習データセット内の最も近い 50 個のアンサンブルメンバーを見つけ、50 個のアンサンブルメンバーに基づいて確率分布を当てはめ、分布の中央値を RUL の推定値として使用します。

mdl = residualSimilarityModel(...
    'Method', 'poly2',...
    'Distance', 'absolute',...
    'NumNearestNeighbors', 50,...
    'Standardize', 1);

fit(mdl, trainDataFused);

性能の評価

類似性 RUL モデルを評価するには、50%、70%、90% のサンプル検証データを使用してその RUL を予測します。

breakpoint = [0.5, 0.7, 0.9];
validationDataTmp = validationDataFused{3}; % use one validation data for illustration

寿命の 50% である最初のブレークポイントより前の検証データを使用します。

bpidx = 1;
validationDataTmp50 = validationDataTmp(1:ceil(end*breakpoint(bpidx)),:);
trueRUL = length(validationDataTmp) - length(validationDataTmp50);
[estRUL, ciRUL, pdfRUL] = predictRUL(mdl, validationDataTmp50);

50% で打ち切られた検証データとその最近傍を可視化します。

figure
compare(mdl, validationDataTmp50);

真の RUL と比較した推定 RUL と、推定 RUL の確率分布を可視化します。

figure
helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

マシンが中間の健全性段階にある場合、推定された RUL と真の RUL の間には比較的大きな誤差があります。この例では、最も類似した 10 個の曲線が最初は近接していますが、故障状態に近づくにつれて分岐し、結果として RUL 分布でおよそ 2 つの最頻値となります。

寿命の 70% である 2 番目のブレークポイントより前の検証データを使用します。

bpidx = 2;
validationDataTmp70 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :);
trueRUL = length(validationDataTmp) - length(validationDataTmp70);
[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp70);

figure
compare(mdl, validationDataTmp70);

figure
helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

観測されるデータが増えるにつれて、RUL の推定は改善されます。

寿命の 90% である 3 番目のブレークポイントより前の検証データを使用します。

bpidx = 3;
validationDataTmp90 = validationDataTmp(1:ceil(end*breakpoint(bpidx)), :);
trueRUL = length(validationDataTmp) - length(validationDataTmp90);
[estRUL,ciRUL,pdfRUL] = predictRUL(mdl, validationDataTmp90);

figure
compare(mdl, validationDataTmp90);

figure
helperPlotRULDistribution(trueRUL, estRUL, pdfRUL, ciRUL)

マシンが故障に近づくと、この例では RUL の推定がさらに改善されます。

ここで、同じ評価手順を検証データセット全体に対して繰り返し、推定された RUL と真の RUL の間の誤差を各ブレークポイントについて計算します。

numValidation = length(validationDataFused);
numBreakpoint = length(breakpoint);
error = zeros(numValidation, numBreakpoint);

for dataIdx = 1:numValidation
    tmpData = validationDataFused{dataIdx};
    for bpidx = 1:numBreakpoint
        tmpDataTest = tmpData(1:ceil(end*breakpoint(bpidx)), :);
        trueRUL = length(tmpData) - length(tmpDataTest);
        [estRUL, ~, ~] = predictRUL(mdl, tmpDataTest);
        error(dataIdx, bpidx) = estRUL - trueRUL;
    end
end

各ブレークポイントの誤差のヒストグラムを、その確率分布とともに可視化します。

[pdf50, x50] = ksdensity(error(:, 1));
[pdf70, x70] = ksdensity(error(:, 2));
[pdf90, x90] = ksdensity(error(:, 3));

figure
ax(1) = subplot(3,1,1);
hold on
histogram(error(:, 1), 'BinWidth', 5, 'Normalization', 'pdf')
plot(x50, pdf50)
hold off
xlabel('Prediction Error')
title('RUL Prediction Error using first 50% of each validation ensemble member')

ax(2) = subplot(3,1,2);
hold on
histogram(error(:, 2), 'BinWidth', 5, 'Normalization', 'pdf')
plot(x70, pdf70)
hold off
xlabel('Prediction Error')
title('RUL Prediction Error using first 70% of each validation ensemble member')

ax(3) = subplot(3,1,3);
hold on
histogram(error(:, 3), 'BinWidth', 5, 'Normalization', 'pdf')
plot(x90, pdf90)
hold off
xlabel('Prediction Error')
title('RUL Prediction Error using first 90% of each validation ensemble member')
linkaxes(ax)

予測誤差を箱ひげ図にプロットして中央値、25 ～ 75 分位数、および外れ値を可視化します。

figure
boxplot(error, 'Labels', {'50%', '70%', '90%'})
ylabel('Prediction Error')
title('Prediction error using different percentages of each validation ensemble member')

予測誤差の平均と標準偏差を計算して可視化します。

errorMean = mean(error)

errorMean = 1×3

   -1.1217    9.5186    9.6540

errorMedian = median(error)

errorMedian = 1×3

    1.3798   11.8172   10.3457

errorSD = std(error)

errorSD = 1×3

   21.7315   13.5194   12.3083

figure
errorbar([50 70 90], errorMean, errorSD, '-o', 'MarkerEdgeColor','r')
xlim([40, 100])
xlabel('Percentage of validation data used for RUL prediction')
ylabel('Prediction Error')
legend('Mean Prediction Error with 1 Standard Deviation Eror bar', 'Location', 'south')

観察されるデータが増えるにつれて、誤差は 0 付近に集中する (外れ値が少ない) ことが示されています。

参考文献

[1] A. Saxena and K. Goebel (2008)."PHM08 Challenge Data Set", NASA Ames Prognostics Data Repository (http://ti.arc.nasa.gov/project/prognostic-data-repository), NASA Ames Research Center, Moffett Field, CA

[2] Roemer, Michael J., Gregory J. Kacprzynski, and Michael H. Schoeller."Improved diagnostic and prognostic assessments using health management information fusion."AUTOTESTCON Proceedings, 2001.IEEE Systems Readiness Technology Conference.IEEE, 2001.

[3] Goebel, Kai, and Piero Bonissone."Prognostic information fusion for constant load systems."Information Fusion, 2005 8th International Conference on.Vol. 2.IEEE, 2005.

[4] Wang, Peng, and David W. Coit."Reliability prediction based on degradation modeling for systems with multiple degradation measures."Reliability and Maintainability, 2004 Annual Symposium-RAMS.IEEE, 2004.

[5] Jardine, Andrew KS, Daming Lin, and Dragan Banjevic."A review on machinery diagnostics and prognostics implementing condition-based maintenance."Mechanical systems and signal processing 20.7 (2006):1483-1510.

補助関数

function data = regimeNormalization(data, centers, centerstats)
% Perform normalization for each observation (row) of the data
% according to the cluster the observation belongs to.
conditionIdx = 3:5;
dataIdx = 6:26;

% Perform row-wise operation
data{:, dataIdx} = table2array(...
    rowfun(@(row) localNormalize(row, conditionIdx, dataIdx, centers, centerstats), ...
    data, 'SeparateInputs', false));
end

function rowNormalized = localNormalize(row, conditionIdx, dataIdx, centers, centerstats)
% Normalization operation for each row.

% Get the operating points and sensor measurements
ops = row(1, conditionIdx);
sensor = row(1, dataIdx);

% Find which cluster center is the closest
dist = sum((centers - ops).^2, 2);
[~, idx] = min(dist);

% Normalize the sensor measurements by the mean and standard deviation of the cluster.
% Reassign NaN and Inf to 0.
rowNormalized = (sensor - centerstats.Mean{idx, :}) ./ centerstats.SD{idx, :};
rowNormalized(isnan(rowNormalized) | isinf(rowNormalized)) = 0;
end

function dataFused = degradationSensorFusion(data, sensorToFuse, weights)
% Combine measurements from different sensors according 
% to the weights, smooth the fused data and offset the data
% so that all the data start from 1

% Fuse the data according to weights
dataToFuse = data{:, cellstr(sensorToFuse)};
dataFusedRaw = dataToFuse*weights;

% Smooth the fused data with moving mean
stepBackward = 10;
stepForward = 10;
dataFused = movmean(dataFusedRaw, [stepBackward stepForward]);

% Offset the data to 1
dataFused = dataFused + 1 - dataFused(1);
end

参考

residualSimilarityModel