Набор данных

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

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

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}$та функция ${\mathit{x}}_{\mathit{i}}$ вычисляется как

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

Размерность и слияние функций

Анализ основных компонентов (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$является lognormal-распределенным и $\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] Bechhoefer, Eric, Brandon Van Hecke, and David He. «Обработка для улучшенного спектрального анализа». Ежегодная конференция Общества прогностики и управления здоровьем, Новый Орлеан, LA, октябрь 2013 года.

[2] Ali, Jaouher Ben, et al. «Оперативная автоматическая диагностика прогрессирующих деградаций подшипников ветряного двигателя в реальных экспериментальных условиях на основе неподконтрольного машинного обучения». Прикладная акустика 132 (2018): 167-181.

[3] Saidi, Lotfi, et al. «Прогноз состояния высокоскоростных подшипников вала ветряного двигателя через спектральные индексы, полученные из Куртоза, и SVR». Прикладная акустика 120 (2017): 1-8.

[4] Кобл, Джейми Баалис. «Объединение источников данных для предсказания оставшегося срока полезного использования - автоматизированный метод для идентификации прогностических параметров». (2010).

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