ドキュメンテーション

このページは前リリースの情報です。該当の英語のページはこのリリースで削除されています。

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

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

データの準備

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

degradationData = helperLoadData('train.txt');
degradationData(1:5)
ans = 5×1 cell array
    {223×26 table}
    {164×26 table}
    {150×26 table}
    {159×26 table}
    {357×26 table}

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

head(degradationData{1})
ans=8×26 table
    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         10.005          0.2501            20          489.05      604.13      1499.5        1310      10.52       15.49       394.88      2318.9      8770.2       1.26          45.4       372.15       2388.1       8120.8       8.6216        0.03          368         2319          100         28.58       17.174  
    1      2         0.0015          0.0003           100          518.67      642.13      1584.5        1404      14.62       21.61       553.67        2388      9045.8        1.3         47.29       521.81       2388.2       8132.9       8.3907        0.03          391         2388          100         38.99       23.362  
    1      3         34.999          0.8401            60          449.44      555.42      1368.2      1122.5       5.48           8       194.93      2222.9      8343.9       1.02         41.92       183.26       2387.9       8063.8       9.3557        0.02          334         2223          100         14.83       8.8555  
    1      4         20.003          0.7005             0          491.19      607.03      1488.4      1249.2       9.35       13.65       334.82      2323.8      8721.5       1.08         44.26       314.84       2388.1       8052.3       9.2231        0.02          364         2324          100         24.42       14.783  
    1      5         42.004          0.8405            40             445      549.52      1354.5      1124.3       3.91        5.71       138.24      2211.8      8314.6       1.02         41.79       130.44       2387.9       8083.7       9.2986        0.02          330         2212          100         10.99       6.4025  
    1      6         20.003          0.7017             0          491.19      607.37      1480.5      1258.9       9.35       13.65       334.51      2323.9      8711.4       1.08          44.4       315.36       2388.1       8053.2       9.2276        0.02          364         2324          100         24.44       14.702  
    1      7             42            0.84            40             445      549.57      1354.4      1131.4       3.91        5.71       139.11      2211.8      8316.9       1.02         42.09       130.16       2387.9         8082       9.3753        0.02          331         2212          100         10.53       6.4254  
    1      8         0.0011               0           100          518.67      642.08      1589.5      1407.6      14.62       21.61       553.48      2388.1      9050.4        1.3          47.5       521.74         2388       8133.3       8.4339        0.03          391         2388          100         38.98       23.234  

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

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.279547.
Replicate 2, 1 iterations, total sum of distances = 0.279547.
Replicate 3, 1 iterations, total sum of distances = 0.279547.
Replicate 4, 1 iterations, total sum of distances = 0.279547.
Replicate 5, 1 iterations, total sum of distances = 0.279547.
Best total sum of distances = 0.279547

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

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.92      1502.1      1311.4      10.52       15.493      394.32        2319      8784.1         1.26      45.496       371.44       2388.2       8133.9       8.6653          0.03      369.74        2319           100       28.525       17.115  
     518.67      642.71      1590.7      1409.4      14.62        21.61       553.3      2388.1      9062.3          1.3      47.557       521.36       2388.1       8140.9       8.4442          0.03      393.27        2388           100       38.809       23.285  
     462.54      536.87      1262.8      1050.6       7.05       9.0275       175.4      1915.4      8014.9      0.93989      36.808       164.56       2028.3         7878       10.916          0.02      307.39        1915         84.93       14.262       8.5552  
        445      549.72      1354.7      1128.2       3.91       5.7158      138.62        2212        8327       1.0202      42.163       130.53         2388       8088.4       9.3775          0.02      331.12        2212           100       10.584       6.3515  
     491.19      607.59        1486      1253.6       9.35       13.657      334.46        2324      8729.1       1.0777      44.466       314.85       2388.2         8065       9.2347      0.022299      365.45        2324           100       24.447        14.67  
     449.44      555.82      1366.9      1131.9       5.48       8.0003      194.43        2223      8355.2       1.0203      41.995       183.01       2388.1       8071.1       9.3344          0.02      334.29        2223           100       14.827       8.8966  

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.4553e-11    0.47617      5.8555      8.3464     1.1618e-12    0.0047272    0.65536     0.094487     18.057       1.51e-13     0.25089      0.52838     0.096183      16.012      0.037598     9.9929e-16     1.4723          0                 0     0.14522     0.086753 
     4.059e-11    0.48566      5.9258      8.8223     3.7307e-14    0.0011406    0.86236     0.068654     20.061     1.2103e-13     0.25937      0.71239     0.069276      17.212      0.036504     9.9929e-16     1.5046          0                 0     0.17565      0.10515 
    2.7685e-11    0.35468      5.2678      6.9664      2.043e-14    0.0043301    0.45074       0.2743     14.741      0.0010469     0.21539      0.33985      0.28932      13.624      0.044142     8.6744e-16     1.3083          0        5.9833e-12     0.11212      0.06722 
             0    0.44169      5.6853      7.5741     1.9763e-13    0.0049401    0.44198      0.30713     18.389      0.0014951     0.23584       0.3432      0.33144      16.917      0.037135      2.231e-15     1.4174          0                 0     0.10778     0.063991 
    4.4456e-11    0.46992      5.7664      7.8679     8.9892e-13    0.0046633    0.59984      0.13032     17.983      0.0042388      0.2362      0.49055      0.13285      15.792      0.038156      0.0042084     1.4428          0                 0     0.13352     0.079731 
    4.3489e-11    0.44341      5.7224      7.4842     3.7485e-13    0.0017642     0.4734      0.28889     17.608      0.0017978     0.23242      0.38243         0.31       15.91      0.038624     8.5009e-16     1.3974          0                 0     0.11311     0.069348 

各領域の統計値を使用して学習データを正規化することができます。各アンサンブル メンバーについて、各行の操作点を抽出し、その各クラスター中心までの距離を計算して、最も近いクラスター中心を見つけます。次に、各センサー測定値について、平均を減算し、そのクラスターの標準偏差で除算します。標準偏差が 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 に変化します。

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

健全性条件 ~ 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.0308
   -0.0308
   -0.0536
    0.0033
   -0.0639
    0.0051
   -0.0408
   -0.0382

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

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)=||yj-yˆj,i||1

ここで、yj はマシン j の健全性インジケーターであり、yˆj,i はマシン i で同定された 2 次多項式モデルを使用する、マシン j の推定された健全性インジケーターです。

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

score(i,j)=exp(-|d(i,j)|)

検証データセット内の 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

   -5.8944    3.1359    3.3555

errorMedian = median(error)
errorMedian = 1×3

   -4.8538    5.3763    3.6580

errorSD = std(error)
errorSD = 1×3

   26.4916   20.0720   18.0313

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

参考

関連するトピック