Мультиклассовое обнаружение неисправностей с использованием смоделированных данных

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

Настройка модели

Этот пример использует много вспомогательных файлов, которые хранятся в zip-файле. Разархивируйте файл, чтобы получить доступ к вспомогательным файлам, загрузить параметры модели и создать библиотеку поршневого насоса.

if ~exist('+mech_hydro_forcesPS','dir')
    unzip('pdmRecipPump_supportingfiles.zip')
end

% Load Parameters
pdmRecipPump_Parameters %Pump
CAT_Pump_1051_DataFile_imported %CAD

% Create Simscape library if needed
if exist('mech_hydro_forcesPS_Lib','file')~=4
    ssc_build mech_hydro_forcesPS
end

Модель поршневого насоса

Поршневой насос состоит из электродвигателя, корпуса насоса, кривошипа насоса и поршней насоса.

mdl = 'pdmRecipPump';
open_system(mdl)

open_system([mdl,'/Pump'])

Модель насоса сконфигурирована, чтобы смоделировать три типа отказов; цилиндрические утечки, блокированный вход и увеличенное трение подшипника. Эти отказы параметризованы как переменные рабочей области и сконфигурированы через диалоговое окно блока насоса.

Симуляция отказа и здоровых данных

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

% Define fault parameter variations
numParValues = 10;
leak_area_set_factor = linspace(0.00,0.036,numParValues);
leak_area_set = leak_area_set_factor*TRP_Par.Check_Valve.In.Max_Area;
leak_area_set = max(leak_area_set,1e-9); % Leakage area cannot be 0
blockinfactor_set = linspace(0.8,0.53,numParValues);
bearingfactor_set = linspace(0,6e-4,numParValues);

Модель насоса сконфигурирована, чтобы включать шум, таким образом запущение модели с теми же значениями параметров отказа приведет к различной симуляции выходные параметры. Это полезно для разработки классификатора, когда это означает, что может быть несколько результатов симуляции для того же условия отказа и серьезности. Чтобы конфигурировать моделирование для таких результатов, создайте векторы значений параметров отказа, где значения не представляют отказов, одиночного отказа, комбинаций двух отказов и комбинаций трех отказов. Для каждой группы (никакой отказ, одиночный отказ, и т.д.) создают 125 комбинаций значений отказа от значений параметров отказа, заданных выше. Это дает в общей сложности 1 000 комбинаций значений параметров отказа.

nPerGroup = 125; % Number of elements in each fault group
    
% No fault simulations
leakArea = repmat(leak_area_set(1),nPerGroup,1);
blockingFactor = repmat(blockinfactor_set(1),nPerGroup,1);
bearingFactor = repmat(bearingfactor_set(1),nPerGroup,1);

% Single fault simulations
idx = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; leak_area_set(idx)'];
blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)];
bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)];
idx = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)];
blockingFactor = [blockingFactor;blockinfactor_set(idx)'];
bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)];
idx = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)];
blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)];
bearingFactor = [bearingFactor;bearingfactor_set(idx)'];

% Double fault simulations
idxA = ceil(10*rand(nPerGroup,1));
idxB = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; leak_area_set(idxA)'];
blockingFactor = [blockingFactor;blockinfactor_set(idxB)'];
bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)];
idxA = ceil(10*rand(nPerGroup,1));
idxB = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; leak_area_set(idxA)'];
blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)];
bearingFactor = [bearingFactor;bearingfactor_set(idxB)'];
idxA = ceil(10*rand(nPerGroup,1));
idxB = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)];
blockingFactor = [blockingFactor;blockinfactor_set(idxA)'];
bearingFactor = [bearingFactor;bearingfactor_set(idxB)'];

% Triple fault simulations
idxA = ceil(10*rand(nPerGroup,1));
idxB = ceil(10*rand(nPerGroup,1));
idxC = ceil(10*rand(nPerGroup,1));
leakArea = [leakArea; leak_area_set(idxA)'];
blockingFactor = [blockingFactor;blockinfactor_set(idxB)'];
bearingFactor = [bearingFactor;bearingfactor_set(idxC)'];

Используйте комбинации параметра отказа, чтобы создать Simulink.SimulationInput объекты. Для каждой симуляции вход гарантирует, что случайный seed собирается по-другому сгенерировать различные результаты.

for ct = numel(leakArea):-1:1
    simInput(ct) = Simulink.SimulationInput(mdl);
    simInput(ct) = setVariable(simInput(ct),'leak_cyl_area_WKSP',leakArea(ct));
    simInput(ct) = setVariable(simInput(ct),'block_in_factor_WKSP',blockingFactor(ct));
    simInput(ct) = setVariable(simInput(ct),'bearing_fault_frict_WKSP',bearingFactor(ct));
    simInput(ct) = setVariable(simInput(ct),'noise_seed_offset_WKSP',ct-1);
end

Используйте generateSimulationEnsemble функционируйте, чтобы запустить симуляции, заданные Simulink.SimulationInput объекты, заданные выше и хранилище результаты в локальной подпапке. Затем создайте simulationEnsembleDatastore от сохраненных результатов.

Обратите внимание на то, что выполнение этих 1 000 параллельных симуляций занимает приблизительно час на стандартном рабочем столе и генерирует приблизительно 620 МБ данных. Возможность, чтобы только запустить первые 10 симуляций предоставляется для удобства.

% Run the simulation and create an ensemble to manage the simulation results
runAll = true;
if runAll
    [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);
else
    [ok,e] = generateSimulationEnsemble(simInput(1:10),fullfile('.','Data')); %#ok<UNRCH>
end
[09-Apr-2018 09:01:38] Checking for availability of parallel pool...
[09-Apr-2018 09:01:38] Loading Simulink on parallel workers...
Analyzing and transferring files to the workers ...done.
[09-Apr-2018 09:01:38] Configuring simulation cache folder on parallel workers...
[09-Apr-2018 09:01:38] Running SetupFcn on parallel workers...
[09-Apr-2018 09:01:39] Loading model on parallel workers...
[09-Apr-2018 09:01:39] Transferring base workspace variables used in the model to parallel workers...
[09-Apr-2018 09:01:41] Running simulations...
[09-Apr-2018 09:02:28] Completed 1 of 1000 simulation runs
[09-Apr-2018 09:02:33] Completed 2 of 1000 simulation runs
[09-Apr-2018 09:02:37] Completed 3 of 1000 simulation runs
[09-Apr-2018 09:02:41] Completed 4 of 1000 simulation runs
[09-Apr-2018 09:02:46] Completed 5 of 1000 simulation runs
[09-Apr-2018 09:02:49] Completed 6 of 1000 simulation runs
[09-Apr-2018 09:02:54] Completed 7 of 1000 simulation runs
[09-Apr-2018 09:02:58] Completed 8 of 1000 simulation runs
[09-Apr-2018 09:03:01] Completed 9 of 1000 simulation runs...
ens = simulationEnsembleDatastore(fullfile('.','Data'));

Обработка и извлечение функций от результатов симуляции

Модель сконфигурирована, чтобы регистрировать давление выхода насоса, вывести поток, частоту вращения двигателя и моторный ток.

ens.DataVariables
ans = 8×1 string array
    "SimulationInput"
    "SimulationMetadata"
    "iMotor_meas"
    "pIn_meas"
    "pOut_meas"
    "qIn_meas"
    "qOut_meas"
    "wMotor_meas"

Поскольку каждый участник в ансамбле предварительно обрабатывает выход насоса, текут и вычисляют индикаторы состояния на основе потока выхода насоса. Индикаторы состояния позже используются в классификации отказов. Поскольку предварительная обработка удаляет первые 0,8 секунды выходного потока, когда это содержит переходные процессы от запуска насоса и симуляции. Как часть предварительной обработки вычисляют спектр мощности выходного потока и используют SimulationInput, чтобы возвратить значения переменных отказа.

Сконфигурируйте ансамбль так, чтобы возвраты только для чтения переменные интереса и вызвали preprocess функция, которая задана в конце этого примера.

ens.SelectedVariables = ["qOut_meas", "SimulationInput"];
reset(ens)
data = read(ens)
data=1×2 table
        qOut_meas                SimulationInput        
    __________________    ______________________________

    [2001×1 timetable]    [1×1 Simulink.SimulationInput]

[flow,flowP,flowF,faultValues] = preprocess(data);

Постройте спектр потока и поток. Отображенные на графике данные являются для отказа свободным условием.

% Figure with nominal
subplot(211);
plot(flow.Time,flow.Data);
subplot(212);
semilogx(flowF,pow2db(flowP));
xlabel('Hz')

Спектр потока показывает резонирующий peaks на ожидаемых частотах. А именно, скорость электродвигателя насоса составляет 950 об/мин или 15,833 Гц, и поскольку насос имеет 3 цилиндра, поток, как ожидают, будет иметь основной принцип на уровне 3*15.833 Гц или 47,5 Гц, а также гармоники во множителях 47,5 Гц. Спектр потока ясно показывает ожидаемый резонирующий peaks. Отказы в одном цилиндре насоса приведут к резонансам на скорости электродвигателя насоса, 15,833 Гц и ее гармониках.

Спектр потока и медленный сигнал дают некоторое представление о возможных индикаторах состояния. А именно, общая статистика сигнала, такая как среднее значение, отклонение, и т.д. а также характеристики спектра. Вычисляются индикаторы состояния спектра, относящиеся к ожидаемым гармоникам, таким как частота с пиковой величиной, энергия приблизительно 15,833 Гц, энергия приблизительно 47,5 Гц, энергия выше 100 Гц. Частота спектрального пика эксцесса также вычисляется.

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

ens.DataVariables = [ens.DataVariables; ...
    "fPeak"; "pLow"; "pMid"; "pHigh"; "pKurtosis"; ...
    "qMean"; "qVar"; "qSkewness"; "qKurtosis"; ...
    "qPeak2Peak"; "qCrest"; "qRMS"; "qMAD"; "qCSRange"];
ens.ConditionVariables = ["LeakFault","BlockingFault","BearingFault"];

feat = extractCI(flow,flowP,flowF);
dataToWrite = [faultValues, feat];
writeToLastMemberRead(ens,dataToWrite{:})

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

%Figure with nominal and faults
figure,
subplot(211);
lFlow = plot(flow.Time,flow.Data,'LineWidth',2);
subplot(212);
lFlowP = semilogx(flowF,pow2db(flowP),'LineWidth',2);
xlabel('Hz')
ct = 1;
lColors = get(lFlow.Parent,'ColorOrder');
lIdx = 2;

% Loop over all members in the ensemble, preprocess 
% and compute the features for each member
while hasdata(ens)
    
    % Read member data
    data = read(ens);
    
    % Preprocess and extract features from the member data
    [flow,flowP,flowF,faultValues] = preprocess(data);
    feat = extractCI(flow,flowP,flowF);
    
    % Add the extracted feature values to the member data
    dataToWrite = [faultValues, feat];
    writeToLastMemberRead(ens,dataToWrite{:})
    
    % Plot member signal and spectrum for every 100th member
    if mod(ct,100) == 0
        line('Parent',lFlow.Parent,'XData',flow.Time,'YData',flow.Data, ...
            'Color', lColors(lIdx,:));
        line('Parent',lFlowP.Parent,'XData',flowF,'YData',pow2db(flowP), ...
            'Color', lColors(lIdx,:));
        if lIdx == size(lColors,1)
            lIdx = 1;
        else
            lIdx = lIdx+1;
        end
    end
    ct = ct + 1;
end

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

Обнаружьте и классифицируйте отказы насоса

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

Сконфигурируйте ансамбль симуляции, чтобы считать индикаторы состояния, и использовать высокое и собрать команды, чтобы загрузить все индикаторы состояния и дать сбой значения переменных в память

% Get data to design a classifier.
reset(ens)
ens.SelectedVariables = [...
    "fPeak","pLow","pMid","pHigh","pKurtosis",...
    "qMean","qVar","qSkewness","qKurtosis",...
    "qPeak2Peak","qCrest","qRMS","qMAD","qCSRange",...
    "LeakFault","BlockingFault","BearingFault"];
idxLastFeature = 14;

% Load the condition indicator data into memory
data = gather(tall(ens));
Starting parallel pool (parpool) using the 'local' profile ...
Preserving jobs with IDs: 1 2 because they contain crash dump files.
You can use 'delete(myCluster.Jobs)' to remove all jobs created with profile local. To create 'myCluster' use 'myCluster = parcluster('local')'.
connected to 6 workers.
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 45 sec
Evaluation completed in 45 sec
data(1:10,:)
ans=10×17 table
    fPeak      pLow       pMid     pHigh     pKurtosis    qMean      qVar     qSkewness    qKurtosis    qPeak2Peak    qCrest     qRMS      qMAD     qCSRange    LeakFault    BlockingFault    BearingFault
    ______    _______    ______    ______    _________    ______    ______    _________    _________    __________    ______    ______    ______    ________    _________    _____________    ____________

    43.909    0.86472    117.63    18.874     276.49      35.572    7.5242    -0.72832      2.7738        13.835      1.1494    35.677    2.2326     42690         1e-09          0.8              0      
    43.909    0.44477    125.92    18.899     12.417      35.576     7.869     -0.7094      2.6338        13.335      1.1449    35.686    2.3204     42697         1e-09          0.8              0      
    43.909     1.1782    137.99    17.526     11.589      35.573    7.4367    -0.72208      2.7136        12.641      1.1395    35.678    2.2407     42695         1e-09          0.8              0      
    44.151     156.74    173.84    21.073      199.5      33.768    12.466    -0.30256      2.4782        17.446      1.2138    33.952    2.8582     40518       2.4e-06          0.8              0      
    43.848    0.71756    110.92    18.579     197.02      35.563    7.5781    -0.72377       2.793         14.14      1.1504    35.669    2.2671     42682         1e-09          0.8              0      
    43.909    0.43673    119.56    20.003     11.589       35.57    7.5028    -0.74797      2.7913        13.833      1.1551    35.676    2.2442     42689         1e-09          0.8              0      
    43.788    0.31617     135.3    19.724     476.82      35.568    7.4406    -0.70964      2.6884        14.685      1.1473    35.673    2.2392     42687         1e-09          0.8              0      
    43.848    0.72747    121.63    19.733     11.589      35.523     7.791    -0.72736      2.7864        14.043      1.1469    35.633    2.2722     42633         1e-09          0.8              0      
    43.848    0.62777    128.85    19.244     11.589      35.541    7.5698     -0.6953      2.6942        13.451      1.1415    35.647    2.2603     42654         1e-09          0.8              0      
    43.848     0.4631    134.83    18.918     12.417      35.561    7.8607    -0.68417      2.6664        13.885      1.1504    35.671    2.3078     42681         1e-09          0.8              0      

Значения переменных отказа для каждого члена ансамбля (строка в таблице данных) могут быть преобразованы во флаги отказа и флаги отказа, объединенные к одному флагу, который получает различное состояние отказа каждого участника.

% Convert the fault variable values to flags
data.LeakFlag = data.LeakFault > 1e-6;
data.BlockingFlag = data.BlockingFault < 0.8;
data.BearingFlag = data.BearingFault > 0; 
data.CombinedFlag = data.LeakFlag+2*data.BlockingFlag+4*data.BearingFlag;

Создайте классификатор, который берет в качестве входа индикаторы состояния и возвращает объединенный флаг отказа. Обучите машину опорных векторов, которая использует 2-е ядро полинома порядка. Используйте cvpartition команда, чтобы разделить членов ансамбля в набор для обучения и набор для валидации.

rng('default') % for reproducibility
predictors = data(:,1:idxLastFeature); 
response = data.CombinedFlag;
cvp = cvpartition(size(predictors,1),'KFold',5);

% Create and train the classifier
template = templateSVM(...
    'KernelFunction', 'polynomial', ...
    'PolynomialOrder', 2, ...
    'KernelScale', 'auto', ...
    'BoxConstraint', 1, ...
    'Standardize', true);
combinedClassifier = fitcecoc(...
    predictors(cvp.training(1),:), ...
    response(cvp.training(1),:), ...
    'Learners', template, ...
    'Coding', 'onevsone', ...
    'ClassNames', [0; 1; 2; 3; 4; 5; 6; 7]);

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

% Check performance by computing and plotting the confusion matrix
actualValue = response(cvp.test(1),:);
predictedValue = predict(combinedClassifier, predictors(cvp.test(1),:));
confdata = confusionmat(actualValue,predictedValue);
figure,
labels = {'None', 'Leak','Blocking', 'Leak & Blocking', 'Bearing', ...
    'Bearing & Leak', 'Bearing & Blocking', 'All'};
h = heatmap(confdata, ...
    'YLabel', 'Actual leak fault', ...
    'YDisplayLabels', labels, ...
    'XLabel', 'Predicted fault', ...
    'XDisplayLabels', labels, ...
    'ColorbarVisible','off');

График беспорядка показывает для каждой комбинации отказов число раз, комбинация отказа была правильно предсказана (диагональные элементы графика) и число раз, комбинация отказа была неправильно предсказана (недиагональные записи).

График беспорядка показывает, что классификатор правильно не классифицировал некоторые условия отказа (от диагональных терминов). Однако никакое условие отказа не было правильно предсказано. В нескольких местах не было предсказано никакое условие отказа, когда был отказ (первый столбец), в противном случае отказ был предсказан несмотря на то, что это не может быть точно правильное условие отказа. В целом точность валидации составляла 84% и точность при предсказании, что существует отказ 98%.

% Compute the overall accuracy of the classifier
sum(diag(confdata))/sum(confdata(:))
ans = 0.6150
% Compute the accuracy of the classifier at predicting 
% that there is a fault
1-sum(confdata(2:end,1))/sum(confdata(:))
ans = 0.9450

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

vData = data(cvp.test(1),:);
b1 = (actualValue==2) & (predictedValue==0);
fData = vData(b1,15:17)
fData=11×3 table
    LeakFault    BlockingFault    BearingFault
    _________    _____________    ____________

      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.71              0      
      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.62              0      
      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.71              0      
      8e-07          0.74              0      
      1e-09          0.74              0      

Найдите случаи в данных о валидации, где фактический отказ был отказом подшипника, но никакой отказ не был предсказан.

b2 = (actualValue==4) & (predictedValue==0);
vData(b2,15:17)
ans =

  0×3 empty table

При исследовании случаев, где никакой отказ не был прогнозирующим, но действительно существовал отказ, показывает, что они происходят, когда блокирующееся значение отказа 0,77 близко к его номинальной стоимости 0,8, или значение отказа подшипника 6.6e-5 близко к его номинальной стоимости 0. Графический вывод спектра для случая с маленьким значением отказа блокирования и сравнения с отказом, свободное условие показывает, что спектры являются очень похожим трудным обнаружением создания. Переквалификация классификатора, но включая блокирующееся значение 0,77 как не дает сбой, условие значительно улучшало бы производительность детектора отказа. В качестве альтернативы использование дополнительных измерений насоса могло предоставить больше информации и улучшить способность обнаружить маленькие отказы блокирования.

% Configure the ensemble to only read the flow and fault variable values
ens.SelectedVariables = ["qOut_meas","LeakFault","BlockingFault","BearingFault"];
reset(ens)

% Load the ensemble member data into memory
data = gather(tall(ens));
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 38 sec
Evaluation completed in 38 sec
% Look for the member that was incorrectly predicted, and 
% compute its power spectrum
idx = ...
    data.BlockingFault == fData.BlockingFault(1) & ...
    data.LeakFault == fData.LeakFault(1) & ...
    data.BearingFault == fData.BearingFault(1);
flow1 = data(idx,1);
flow1 = flow1.qOut_meas{1};
[flow1P,flow1F] = pspectrum(flow1);

% Look for a member that does not have any faults
idx = ...
    data.BlockingFault == 0.8 & ...
    data.LeakFault == 1e-9 & ...
    data.BearingFault == 0;
flow2 = data(idx,1);
flow2 = flow2.qOut_meas{1};
[flow2P,flow2F] = pspectrum(flow2);

% Plot the power spectra
semilogx(...
    flow1F,pow2db(flow1P),...
    flow2F,pow2db(flow2P));
xlabel('Hz')
legend('Small blocking fault','No fault')

Заключение

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

Вспомогательные Функции

function [flow,flowSpectrum,flowFrequencies,faultValues] = preprocess(data)
% Helper function to preprocess the logged reciprocating pump data.

% Remove the 1st 0.8 seconds of the flow signal
tMin = seconds(0.8);
flow = data.qOut_meas{1};
flow = flow(flow.Time >= tMin,:);
flow.Time = flow.Time - flow.Time(1);

% Ensure the flow is sampled at a uniform sample rate
flow = retime(flow,'regular','linear','TimeStep',seconds(1e-3));

% Remove the mean from the flow and compute the flow spectrum
fA = flow;
fA.Data = fA.Data - mean(fA.Data);
[flowSpectrum,flowFrequencies] = pspectrum(fA,'FrequencyLimits',[2 250]);

% Find the values of the fault variables from the SimulationInput
simin = data.SimulationInput{1};
vars = {simin.Variables.Name};
idx = strcmp(vars,'leak_cyl_area_WKSP');
LeakFault = simin.Variables(idx).Value;
idx = strcmp(vars,'block_in_factor_WKSP');
BlockingFault = simin.Variables(idx).Value;
idx = strcmp(vars,'bearing_fault_frict_WKSP');
BearingFault = simin.Variables(idx).Value;

% Collect the fault values in a cell array
faultValues = {...
    'LeakFault', LeakFault, ...
    'BlockingFault', BlockingFault, ...
    'BearingFault', BearingFault};
end

function ci = extractCI(flow,flowP,flowF)
% Helper function to extract condition indicators from the flow signal 
% and spectrum.

% Find the frequency of the peak magnitude in the power spectrum.
pMax = max(flowP);
fPeak = flowF(flowP==pMax);

% Compute the power in the low frequency range 10-20 Hz.
fRange = flowF >= 10 & flowF <= 20;
pLow = sum(flowP(fRange));

% Compute the power in the mid frequency range 40-60 Hz.
fRange = flowF >= 40 & flowF <= 60;
pMid = sum(flowP(fRange));

% Compute the power in the high frequency range >100 Hz.
fRange = flowF >= 100;
pHigh = sum(flowP(fRange));

% Find the frequency of the spectral kurtosis peak
[pKur,fKur] = pkurtosis(flow);
pKur = fKur(pKur == max(pKur));

% Compute the flow cumulative sum range.
csFlow = cumsum(flow.Data);
csFlowRange = max(csFlow)-min(csFlow);

% Collect the feature and feature values in a cell array. 
% Add flow statistic (mean, variance, etc.) and common signal 
% characteristics (rms, peak2peak, etc.) to the cell array.
ci = {...
    'qMean', mean(flow.Data), ...
    'qVar',  var(flow.Data), ...
    'qSkewness', skewness(flow.Data), ...
    'qKurtosis', kurtosis(flow.Data), ...
    'qPeak2Peak', peak2peak(flow.Data), ...
    'qCrest', peak2rms(flow.Data), ...
    'qRMS', rms(flow.Data), ...
    'qMAD', mad(flow.Data), ...
    'qCSRange',csFlowRange, ...
    'fPeak', fPeak, ...
    'pLow', pLow, ...
    'pMid', pMid, ...
    'pHigh', pHigh, ...
    'pKurtosis', pKur(1)};
end 

Смотрите также

Похожие темы

Для просмотра документации необходимо авторизоваться на сайте