Набор данных

Набор данных собран из ветряного двигателя на 2 мВт высокоскоростной вал, управляемый механизмом шестерни с 20 зубами [1]. Сигнал вибрации 6 секунд получался каждый день в течение 50 дней подряд (17 марта существует 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-файл и сохраните его в той же директории как этот live скрипт. Разархивируйте файл с помощью этой команды. Временной шаг измерения для полного набора данных составляет 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]
             :                      :                   :
             :                      :                   :

Частота дискретизации сигнала вибрации составляет 97 656 Гц.

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 кГц постепенно увеличивают как условие машины, ухудшается. Статистическими функциями спектрального эксцесса, такими как среднее значение, стандартное отклонение и т.д., будут потенциальные индикаторы ухудшения подшипника [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"];

Поскольку таблица функции мала достаточно, чтобы уместиться в памяти (50 15), собрать таблицу перед обработкой. Для больших данных рекомендуется выполнить операции в высоком формате, пока вы не уверены, что выход мал достаточно, чтобы уместиться в памяти.

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], используется, чтобы определить количество заслуги функций цели прогноза.

Монотонность $\mathit{i}$функция th ${\mathit{x}}_{\mathit{i}}$ вычисляется как

где $\mathit{n}$ количество точек измерения, в этом случае $\mathit{n}=50$. $\mathit{m}$ количество проверенных машин, в этом случае $\mathit{m}=1$. ${\mathit{x}}_{\mathit{i}}^{\mathit{j}}$ $\mathit{i}$функция th, измеренная на $\mathit{j}$машина th. $\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.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
      ⋮

Сокращение размерности и Fusion функции

Анализ главных компонентов (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)

Экспоненциальная модель ухудшения задана как

где $\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);

Если исторические данные доступны, используйте fit метод предоставлен exponentialDegradationModel оценить уголовное прошлое и прерывание. Однако исторические данные не доступны для этого ветряного двигателя, переносящего набор данных. Предшествующие из наклонных параметров выбраны произвольно с большими отклонениями ($\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 также.

Отношение между изменением медицинского индикатора и изменением шума может быть выведено как

Здесь стандартное отклонение шума принято, чтобы вызвать 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 между $\alpha$ связанный истинного RUL вычисляется как показатель производительности модели:

где $\mathit{r}\left(\mathit{t}\right)$ предполагаемый RUL во время $\mathit{t}$, ${\mathit{r}}^{*}\left(\mathit{t}\right)$ истинный RUL во время $\mathit{t}$, $\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] <http://data-acoustics.com/measurements/bearing-faults/bearing-3/> Bechhoefer, Эрик, Брэндон Ван Хек и Дэвид Хэ. "Обрабатывая для улучшенного спектрального анализа". Ежегодная конференция Предзнаменований и медицинского Общества управления, Нового Орлеана, LA, октябрь 2013.

[2] Али, Джэоуэр Бен, и др. "Онлайновый автоматический диагноз подшипников ветряного двигателя прогрессивные ухудшения при действительных экспериментальных условиях на основе безнадзорного машинного обучения". Прикладная Акустика 132 (2018): 167-181.

[3] Саиди, Lotfi, и др. "Ветряной двигатель высокоскоростной медицинский прогноз подшипников вала через спектральные Выведенные из эксцесса индексы и SVR". Прикладная Акустика 120 (2017): 1-8.

[4] Плоскодонная рыбачья лодка, Джейми Баалис. "Объединяя источники данных, чтобы предсказать остающийся срок полезного использования – автоматизированный метод, чтобы идентифицировать предвещающие параметры". (2010).

[5] Saxena, Abhinav, и др. "Метрики для оффлайновой оценки предвещающей производительности". Международный журнал Предзнаменований и медицинского управления 1.1 (2010): 4-23.