풍력 터빈 고속 베어링 예지진단

라이브 스크립트 열기

이 예제에서는 지수 성능 저하 모델을 구축하여 풍력 터빈 베어링의 잔여 수명(RUL)을 실시간으로 예측하는 방법을 보여줍니다. 지수 성능 저하 모델은 파라미터 사전 값과 최신 측정값을 기반으로 RUL을 예측합니다(과거 RTF(run-to-failure) 데이터는 모델의 파라미터 사전 값을 추정하는 데 도움이 되긴 하지만 반드시 필요하지는 않음). 모델은 유의미한 성능 저하 추세를 실시간으로 감지하고 새로운 관측값이 생기면 파라미터 사전 값을 업데이트합니다. 이 예제에서는 데이터 가져오기 및 탐색, 특징 추출 및 전처리, 특징 중요도 순위 지정 및 융합, 모델 피팅 및 예측, 성능 분석이라는 일반적인 예지진단 워크플로를 따릅니다.

데이터셋

데이터셋은 20개 톱니 피니언 기어로 구동되는 2MW 풍력 터빈 고속 축에서 수집되었습니다 [1]. 연속 50일 동안 매일 6초 분량의 진동 신호가 수집되었습니다(3월 17일에는 2건의 측정값이 있는데, 이 예제에서 이는 2일로 취급됨). 50일 기간 안에 내륜 결함이 발생하여 그 결과 베어링 고장이 발생했습니다.

툴박스에서 데이터셋의 간소 버전을 사용할 수 있습니다. 간소 데이터셋을 사용하려면 데이터셋을 현재 폴더에 복사하고 쓰기 권한을 활성화하십시오.

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

간소 데이터셋의 측정 시간 스텝은 5일입니다.

timeUnit = '\times 5 day';

전체 데이터셋을 사용하려면 https://github.com/mathworks/WindTurbineHighSpeedBearingPrognosis-Data에서 전체 리포지토리를 zip 파일로 다운로드하고 이 라이브 스크립트와 동일한 디렉터리에 저장하십시오. 다음 명령을 사용하여 파일의 압축을 풉니다. 전체 데이터셋의 측정 시간 스텝은 1일입니다.

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

이 예제의 결과는 전체 데이터셋에서 생성한 것입니다. 전체 데이터셋을 다운로드하여 이 예제를 실행하는 것이 좋습니다. 간소 데이터셋에서 생성된 결과는 유의미하지 않을 수 있습니다.

데이터 가져오기

풍력 터빈 데이터의 fileEnsembleDatastore를 만듭니다. 데이터는 진동 신호와 회전속도계 신호를 포함합니다. fileEnsembleDatastore는 파일 이름을 구문 분석하고 날짜 정보를 IndependentVariables로 추출합니다. 자세한 내용은 이 예제와 연결된 지원 파일의 헬퍼 함수를 참조하십시오.

hsbearing = fileEnsembleDatastore(...
    fullfile('.', 'WindTurbineHighSpeedBearingPrognosis-Data-master'), ...
    '.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]
             :                      :                   :
             :                      :                   :

진동 신호의 샘플 레이트는 97656Hz입니다.

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 스케일로 정규화된 것입니다. 기계 상태의 성능이 저하됨에 따라 10kHz 부근의 스펙트럼 첨도 값이 점진적으로 증가하는 것을 볼 수 있습니다. 평균, 표준편차 등과 같은 통계 특징은 베어링 성능 저하의 잠재적인 지표가 됩니다 [3].

특징 추출

이전 섹션의 분석을 기반으로 하여 시간 영역 신호에서 도출된 통계 특징들과 스펙트럼 첨도를 추출하려고 합니다. 특징에 대한 자세한 수학적 내용은 [2-3]을 참조하십시오.

먼저 DataVariables에 특징 이름을 미리 할당한 후에 fileEnsembleDatastore에 씁니다.

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"];

특징 테이블이 메모리에 담을 수 있을 정도로 작으므로(50x15) 처리에 앞서 테이블을 수집합니다. 빅데이터의 경우에는 출력값이 메모리에 담을 수 있을 정도로 작다고 확신할 수 있을 때까지 연산을 tall형 형식으로 수행하는 것이 좋습니다.

featureTable = gather(tall(hsbearing));

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

시간 정보가 항상 특징 값과 연결되도록 테이블을 타임테이블로 변환합니다.

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.001253    0.025674     -0.22882        3.362  
    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.0046823    0.020888     0.057651       3.3508  
    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.0011084    0.022705     -0.50004       4.9953  
    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.0087035    0.034456       1.4705       8.1235  
    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.013559    0.032193      0.11355        3.848  
    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.0017539    0.028326      -0.1288       3.8072  
    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.0039353    0.029292      -1.4734       8.1242  
    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.0013107    0.022468      0.27438       2.8597  
    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.0023769    0.047767      -2.5832       20.171  
    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.0076327    0.030418      0.52322       4.0082  
    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.0084907    0.037876       2.3753       11.491  
    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.019624    0.055537       3.1986       17.796  
    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.016315    0.064516       2.8735       11.632  
    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.0011097    0.051539    -0.056774       7.0712  
    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.0040392    0.066254     -0.39587       12.111  
    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.020676    0.077728       2.6038       11.088  
      ⋮

특징 후처리

추출된 특징은 보통 잡음과 연관됩니다. 반대의 추세를 갖는 잡음은 때로는 RUL 예측에 해로울 수 있습니다. 이에 더해 조금 뒤에 소개할 특징 성능 메트릭인 단조 특성은 잡음에 강인하지 않습니다. 따라서 추출된 특징에 지연 윈도우가 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]에서 제시하는 단조 특성을 사용하여 예지진단을 위해 특징의 이점을 수량화합니다.

$i$ 번째 특징 $x_{i}$ 의 단조 특성은 다음과 같이 계산됩니다.

$Monotonicity (x_{i}) = \frac{1}{m} \sum_{j = 1}^{m} \frac{| number of positive diff (x_{i}^{j}) - number of negative diff (x_{i}^{j}) |}{n - 1}$

여기서 $n$ 은 측정 점의 개수로, 여기서는 $n = 50$ 입니다. $m$ 은 모니터링되는 기계의 개수로, 여기서는 $m = 1$ 입니다. $x_{i}^{j}$ 는 $j$ 번째 기계에서 측정된 $i$ 번째 특징입니다. $diff (x_{i}^{j}) = x_{i}^{j} (t) - x_{i}^{j} (t - 1)$ , 즉 신호 $x_{i}^{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.025674
    08-Mar-2013 02:34:21    0.29507     3.0075        1.254          3.5428       0.023281
    09-Mar-2013 02:33:43    0.26962     3.0125        1.254          3.6541       0.023089
    10-Mar-2013 03:01:02    0.25565     3.0197       1.2544          3.7722       0.025931
    11-Mar-2013 03:00:24    0.24756     3.0247       1.2546          3.7793       0.027183
    12-Mar-2013 06:17:10    0.25519     3.0236       1.2544          3.7479       0.027374
    13-Mar-2013 06:34:04      0.233     3.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 = 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);

처음 두 개의 주성분으로 구성된 공간에서 데이터를 시각화합니다.

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

플롯은 기계가 고장에 가까워지면 첫 번째 주성분이 증가함을 나타냅니다. 따라서 첫 번째 주성분은 융합된 건전성 지표로서 적절해 보입니다.

healthIndicator = PCA1;

건전성 지표를 시각화합니다.

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

잔여 수명(RUL) 추정을 위해 지수 성능 저하 모델 피팅하기

지수 성능 저하 모델은 다음과 같이 정의됩니다.

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

여기서 $h (t)$ 는 건전성 지표로, 시간의 함수입니다. $ϕ$ 는 절편 항으로, 상수로 간주됩니다. $θ$ 와 $β$ 는 모델의 기울기를 결정하는 임의 파라미터입니다. 여기서 $θ$ 는 로그 정규분포이고 $β$ 는 가우스 분포입니다. 각 시간 스텝 $t$ 에서 $θ$ 와 $β$ 의 분포는 $h (t)$ 의 최신 관측값을 기준으로 사후 확률로 업데이트됩니다. $ϵ$ 은 $N (0, σ^{2})$ 을 따르는 가우스 백색 잡음입니다. 지수의 $- \frac{σ^{2}}{2}$ 항은 $h (t)$ 의 기대값이 다음을 충족하도록 만들기 위한 것입니다.

$E [h (t) | θ, β] = ϕ + θ \exp (β t)$ .

여기서는 지수 성능 저하 모델이 앞 섹션에서 추출한 건전성 지표에 피팅되고, 다음 섹션에서는 성능이 평가됩니다.

먼저 0에서 시작되도록 건전성 지표를 이동합니다.

healthIndicator = healthIndicator - healthIndicator(1);

임계값은 일반적으로 기계에 대한 과거 기록이나 일부 영역별 지식을 기반으로 선택합니다. 이 데이터셋에 과거 데이터가 없으므로 임계값으로 건전성 지표의 마지막 값을 선택합니다. 평활화로 인한 지연의 영향이 부분적으로 완화될 수 있도록 평활화된 과거 데이터를 기반으로 임계값을 선택하는 것이 좋습니다.

threshold = healthIndicator(end);

과거 데이터가 있는 경우 exponentialDegradationModel에서 제공하는 fit 메서드를 사용하여 사전 값과 절편을 추정합니다. 그러나 이 풍력 터빈 베어링 데이터셋에는 과거 데이터가 없습니다. 기울기 파라미터의 사전 값은 모델이 주로 관측된 데이터에 의존하도록 큰 분산을 사용하여 임의로 선택됩니다( $E (θ) = 1, Var (θ) = 10^{6}, E (β) = 1, Var (β) = 10^{6}$ ). $E [h (0)] = ϕ + E (θ)$ 를 기반으로, 모델도 0에서 시작되도록 절편 $ϕ$ 를 $- 1$ 로 설정합니다.

건전성 지표의 변동과 잡음의 변동 사이의 관계는 다음과 같이 유도할 수 있습니다.

$Δ h (t) \approx (h (t) - ϕ) Δ ϵ (t)$

여기서 잡음의 표준편차는 건전성 지표가 임계값에 가까울 때 건전성 지표 변동의 10%를 유발한다고 가정합니다. 따라서 잡음의 표준편차는 $\frac{10 % \cdot threshold}{threshold - ϕ}$ 로 표현할 수 있습니다.

지수 성능 저하 모델은 기울기의 유의성을 평가하는 기능도 제공합니다. 건전성 지표의 유의미한 기울기가 감지되면 모델은 이전 관측값을 망각하고 원래의 사전 값을 기반으로 추정을 다시 시작합니다. 감지 알고리즘의 민감도는 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 추정의 애니메이션입니다.

성능 분석

$α$ - $λ$ 플롯은 예지진단 성능 분석에 사용되며 [5], 여기서 $α$ 한계는 20%로 설정되었습니다. 추정된 RUL이 실제 RUL의 $α$ 한계 사이에 있을 확률은 모델의 성능 메트릭으로 계산됩니다.

$\Pr (r^{*} (t) - α r^{*} (t) < r (t) < r^{*} (t) + α r^{*} (t) | Θ (t))$

여기서 $r (t)$ 는 시간 $t$ 에서의 추정된 RUL이고, $r^{*} (t)$ 는 시간 $t$ 에서의 실제 RUL이고, $Θ (t)$ 는 시간 $t$ 에서의 추정된 모델 파라미터입니다.

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

사전 설정된 사전 값이 실제 사전 값을 반영하지 않으므로 모델이 올바른 파라미터 분포에 맞게 조정되려면 보통 몇 개의 시간 스텝이 필요합니다. 더 많은 데이터 점이 생김에 따라 예측이 더 정확해집니다.

$α$ 한계 내에서 예측된 RUL의 확률을 시각화합니다.

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.