Набор данных

Набор данных собирается с 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-файла и сохраните его в том же каталоге, что и этот сценарий. Распакуйте файл с помощью этой команды. Шаг времени измерения для полного набора данных составляет 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]
             :                      :                   :
             :                      :                   :

Частота дискретизации сигнала вибрации 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 перед их записью в файл EnsingedDatastore.

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{}}_{\mathit{xi}}$ вычисляется как

где $\mathit{n}$ - количество точек измерения, в данном случае $\mathit{n}\mathrm{=}$50$\mathit{.}$ m - количество контролируемых машин, в данном $\mathit{случае}\mathrm{m}$= ${\mathit{}}_{\mathit{1}}^{\mathit{.}}$ xij - $\mathit{}$i-й признак, измеренный $\mathit{}$на j-ой $\mathrm{машине.}{\mathit{}}_{\mathit{}}^{\mathit{diff}}({\mathit{}}_{\mathit{}}^{\mathit{xij}})\mathit{}={\mathit{}}_{\mathit{}}^{\mathit{xij}}\mathit{(}t\mathrm{)}$-xij (t-1), т.е. разность ${\mathit{}}_{\mathit{}}^{\mathit{}}$сигнала xij.

% 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}\mathit{(}t$) - показатель здоровья как функция времени. $\text{start-}$термин перехвата, рассматриваемый как константа. $$λ и $\text{β}$ - случайные параметры, определяющие наклон модели, где $\lambda$- логнормально-распределённый, а $\beta$- гауссово-распределённый. На каждом временном шаге $\mathit{t}$распределение $\lambda$и $\beta$обновляется до заднего на основе последнего наблюдения $\mathit{h}\mathit{(}t$). $\text{ϵ}$ является гауссовым белым шумом, дающим $\mathit{N}\mathrm{(}0{,}^{\mathrm{}}$start2). $Член\frac{{}^{\mathrm{-}}}{\mathrm{}}$start22 в экспоненциальном является, чтобы сделать ожидание $\mathit{}\mathit{h}$(t) удовлетворяющим

$\mathit{E}\mathit{[}h\mathit{}(t)|\mathrm{}\mathit{}$

Здесь экспоненциальная модель деградации соответствует показателю работоспособности, извлеченному в последнем разделе, и рабочие характеристики оцениваются в следующем разделе.

Сначала переместите индикатор работоспособности таким образом, чтобы он начинался с 0.

healthIndicator = healthIndicator - healthIndicator(1);

Выбор порога обычно основан на исторических записях машины или некоторых знаниях, специфичных для домена. Поскольку в этом наборе данных отсутствуют исторические данные, в качестве порогового значения выбирается последнее значение индикатора работоспособности. Рекомендуется выбирать порог на основе сглаженных (исторических) данных, чтобы эффект задержки сглаживания был частично смягчен.

threshold = healthIndicator(end);

При наличии исторических данных используйте fit способ, предоставляемый exponentialDegradationModel для оценки приоров и перехвата. Однако для данного набора данных по подшипникам ветровых турбин исторические данные отсутствуют. Предшествующие из параметров наклона выбираются произвольно с большими дисперсиями ($\mathit{E}(start)\mathrm{}=\mathrm{1,}Var({\mathrm{start)}}^{\mathrm{}}=\mathit{}106,\mathrm{}E\mathrm{(\beta })={\mathrm{1,}}^{\mathrm{}}$Var (β) = 106), так что модель в основном опирается на наблюдаемые $\mathit{}данные\mathit{.}\mathrm{}\mathit{}На\; основеE$[h (0)] = ϕ + E $$(θ), $\mathrm{}$перехватите ϕ, установлен в-1 так, чтобы модель началась от 0 также.

Взаимосвязь между изменением показателя здоровья и изменением шума может быть получена следующим образом:

Здесь предполагается, что стандартное отклонение шума вызывает 10% изменения индикатора работоспособности, когда он находится вблизи порога. Поэтому стандартное отклонение шума может быть представлено как $\frac{\mathrm{}\mathrm{}}{\mathrm{}\mathrm{10\%\cdot thresholdthreshold-\varphi }\text{.}}$

Экспоненциальная модель ухудшения также обеспечивает функциональные возможности для оценки значимости наклона. Как только обнаружен значительный наклон индикатора работоспособности, модель забудет предыдущие наблюдения и перезапустит оценку на основе исходных предварительных значений. Чувствительность алгоритма обнаружения можно настроить, указав 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], где $\mathrm{\alpha -граница}$установлена в 20%. Вероятность того, что оцененный RUL находится между $\mathrm{\alpha -границей}$ истинного RUL, вычисляется как показатель производительности модели:

где $\mathit{r}\mathit{(}t$) - оцененное значение RUL в момент $\mathit{}$времени t, ${\mathit{}}^r\mathit{*}$(t) - истинное значение RUL в $\mathit{}$момент времени $t\mathit{}$$\mathit{}$,

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

Поскольку предустановленный предшествующий уровень не отражает истинный предшествующий уровень, модель обычно нуждается в нескольких временных шагах, чтобы приспособиться к правильному распределению параметров. Прогнозирование становится более точным по мере появления большего количества точек данных.

Визуализируйте вероятность предсказанного RUL в пределах $\mathrm{\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] Беххёфер, Эрик, Брэндон Ван Хекке и Дэвид Хе. «Обработка для улучшенного спектрального анализа». Ежегодная конференция Общества по прогностике и управлению здравоохранением, Новый Орлеан, Лос-Анджелес, октябрь 2013 года.

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

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

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

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