В этом примере показано, как использовать модель 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 комбинаций значений параметров отказа. Обратите внимание на то, что выполнение этих 1 000 параллельных симуляций занимает приблизительно час на стандартном рабочем столе и генерирует приблизительно 620 МБ данных. Чтобы уменьшать время симуляции, сократите количество комбинаций отказа к 20 путем изменения runAll = true к runAll = false. Note that a больший набор данных приводит к более устойчивому классификатору.
% Set number of elements in each fault group runAll = false; if runAll % Create a large dataset to build a robust classifier nPerGroup = 125; else % Create a smaller dataset to reduce simulation time nPerGroup = 20; %#ok<UNRCH> end % 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 от сохраненных результатов.
% Run the simulation and create an ensemble to manage the simulation % results if isfolder('./Data') % Delete existing mat files delete('./Data/*.mat') end [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);
[08-Sep-2020 21:01:46] Checking for availability of parallel pool... Starting parallel pool (parpool) using the 'local' profile ... Connected to the parallel pool (number of workers: 6). [08-Sep-2020 21:02:27] Starting Simulink on parallel workers... [08-Sep-2020 21:02:55] Configuring simulation cache folder on parallel workers... [08-Sep-2020 21:02:55] Loading model on parallel workers... [08-Sep-2020 21:03:02] Transferring base workspace variables used in the model to parallel workers... [08-Sep-2020 21:03:08] Running simulations... [08-Sep-2020 21:04:21] Completed 1 of 160 simulation runs [08-Sep-2020 21:04:21] Completed 2 of 160 simulation runs [08-Sep-2020 21:04:22] Completed 3 of 160 simulation runs [08-Sep-2020 21:04:23] Completed 4 of 160 simulation runs [08-Sep-2020 21:04:24] Completed 5 of 160 simulation runs [08-Sep-2020 21:04:24] Completed 6 of 160 simulation runs [08-Sep-2020 21:04:47] Completed 7 of 160 simulation runs [08-Sep-2020 21:04:48] Completed 8 of 160 simulation runs [08-Sep-2020 21:04:49] Completed 9 of 160 simulation runs [08-Sep-2020 21:04:50] Completed 10 of 160 simulation runs [08-Sep-2020 21:04:52] Completed 11 of 160 simulation runs [08-Sep-2020 21:04:53] Completed 12 of 160 simulation runs [08-Sep-2020 21:05:14] Completed 13 of 160 simulation runs [08-Sep-2020 21:05:15] Completed 14 of 160 simulation runs [08-Sep-2020 21:05:17] Completed 15 of 160 simulation runs [08-Sep-2020 21:05:18] Completed 16 of 160 simulation runs [08-Sep-2020 21:05:20] Completed 17 of 160 simulation runs [08-Sep-2020 21:05:21] Completed 18 of 160 simulation runs [08-Sep-2020 21:05:41] Completed 19 of 160 simulation runs [08-Sep-2020 21:05:43] Completed 20 of 160 simulation runs [08-Sep-2020 21:05:45] Completed 21 of 160 simulation runs [08-Sep-2020 21:05:47] Completed 22 of 160 simulation runs [08-Sep-2020 21:05:49] Completed 23 of 160 simulation runs [08-Sep-2020 21:05:51] Completed 24 of 160 simulation runs [08-Sep-2020 21:06:08] Completed 25 of 160 simulation runs [08-Sep-2020 21:06:10] Completed 26 of 160 simulation runs [08-Sep-2020 21:06:13] Completed 27 of 160 simulation runs [08-Sep-2020 21:06:16] Completed 28 of 160 simulation runs [08-Sep-2020 21:06:18] Completed 29 of 160 simulation runs [08-Sep-2020 21:06:21] Completed 30 of 160 simulation runs [08-Sep-2020 21:06:34] Completed 31 of 160 simulation runs [08-Sep-2020 21:06:38] Completed 32 of 160 simulation runs [08-Sep-2020 21:06:41] Completed 33 of 160 simulation runs [08-Sep-2020 21:06:44] Completed 34 of 160 simulation runs [08-Sep-2020 21:06:48] Completed 35 of 160 simulation runs [08-Sep-2020 21:06:51] Completed 36 of 160 simulation runs [08-Sep-2020 21:07:01] Completed 37 of 160 simulation runs [08-Sep-2020 21:07:05] Completed 38 of 160 simulation runs [08-Sep-2020 21:07:09] Completed 39 of 160 simulation runs [08-Sep-2020 21:07:13] Completed 40 of 160 simulation runs [08-Sep-2020 21:07:17] Completed 41 of 160 simulation runs [08-Sep-2020 21:07:21] Completed 42 of 160 simulation runs [08-Sep-2020 21:07:32] Completed 43 of 160 simulation runs [08-Sep-2020 21:07:37] Completed 44 of 160 simulation runs [08-Sep-2020 21:07:41] Completed 45 of 160 simulation runs [08-Sep-2020 21:07:46] Completed 46 of 160 simulation runs [08-Sep-2020 21:07:50] Completed 47 of 160 simulation runs [08-Sep-2020 21:07:58] Completed 48 of 160 simulation runs [08-Sep-2020 21:08:04] Completed 49 of 160 simulation runs [08-Sep-2020 21:08:10] Completed 50 of 160 simulation runs [08-Sep-2020 21:08:16] Completed 51 of 160 simulation runs [08-Sep-2020 21:08:23] Completed 52 of 160 simulation runs [08-Sep-2020 21:08:33] Completed 53 of 160 simulation runs [08-Sep-2020 21:08:41] Completed 54 of 160 simulation runs [08-Sep-2020 21:08:51] Completed 55 of 160 simulation runs [08-Sep-2020 21:09:01] Completed 56 of 160 simulation runs [08-Sep-2020 21:09:11] Completed 57 of 160 simulation runs [08-Sep-2020 21:09:22] Completed 58 of 160 simulation runs [08-Sep-2020 21:09:34] Completed 59 of 160 simulation runs [08-Sep-2020 21:09:45] Completed 60 of 160 simulation runs [08-Sep-2020 21:09:57] Completed 61 of 160 simulation runs [08-Sep-2020 21:10:08] Completed 62 of 160 simulation runs [08-Sep-2020 21:10:16] Completed 63 of 160 simulation runs [08-Sep-2020 21:10:25] Completed 64 of 160 simulation runs [08-Sep-2020 21:10:36] Completed 65 of 160 simulation runs [08-Sep-2020 21:10:44] Completed 66 of 160 simulation runs [08-Sep-2020 21:10:57] Completed 67 of 160 simulation runs [08-Sep-2020 21:11:10] Completed 68 of 160 simulation runs [08-Sep-2020 21:11:21] Completed 69 of 160 simulation runs [08-Sep-2020 21:11:31] Completed 70 of 160 simulation runs [08-Sep-2020 21:11:44] Completed 71 of 160 simulation runs [08-Sep-2020 21:11:55] Completed 72 of 160 simulation runs [08-Sep-2020 21:12:11] Completed 73 of 160 simulation runs [08-Sep-2020 21:12:28] Completed 74 of 160 simulation runs [08-Sep-2020 21:12:43] Completed 75 of 160 simulation runs [08-Sep-2020 21:12:59] Completed 76 of 160 simulation runs [08-Sep-2020 21:13:12] Completed 77 of 160 simulation runs [08-Sep-2020 21:13:27] Completed 78 of 160 simulation runs [08-Sep-2020 21:13:43] Completed 79 of 160 simulation runs [08-Sep-2020 21:13:58] Completed 80 of 160 simulation runs [08-Sep-2020 21:14:18] Completed 81 of 160 simulation runs [08-Sep-2020 21:14:30] Completed 82 of 160 simulation runs [08-Sep-2020 21:14:36] Completed 83 of 160 simulation runs [08-Sep-2020 21:14:42] Completed 84 of 160 simulation runs [08-Sep-2020 21:14:47] Completed 85 of 160 simulation runs [08-Sep-2020 21:14:53] Completed 86 of 160 simulation runs [08-Sep-2020 21:14:59] Completed 87 of 160 simulation runs [08-Sep-2020 21:15:05] Completed 88 of 160 simulation runs [08-Sep-2020 21:15:12] Completed 89 of 160 simulation runs [08-Sep-2020 21:15:18] Completed 90 of 160 simulation runs [08-Sep-2020 21:15:26] Completed 91 of 160 simulation runs [08-Sep-2020 21:15:33] Completed 92 of 160 simulation runs [08-Sep-2020 21:15:40] Completed 93 of 160 simulation runs [08-Sep-2020 21:15:46] Completed 94 of 160 simulation runs [08-Sep-2020 21:15:53] Completed 95 of 160 simulation runs [08-Sep-2020 21:16:00] Completed 96 of 160 simulation runs [08-Sep-2020 21:16:07] Completed 97 of 160 simulation runs [08-Sep-2020 21:16:14] Completed 98 of 160 simulation runs [08-Sep-2020 21:16:21] Completed 99 of 160 simulation runs [08-Sep-2020 21:16:29] Completed 100 of 160 simulation runs [08-Sep-2020 21:16:36] Completed 101 of 160 simulation runs [08-Sep-2020 21:16:41] Completed 102 of 160 simulation runs [08-Sep-2020 21:16:48] Completed 103 of 160 simulation runs [08-Sep-2020 21:16:53] Completed 104 of 160 simulation runs [08-Sep-2020 21:16:58] Completed 105 of 160 simulation runs [08-Sep-2020 21:17:02] Completed 106 of 160 simulation runs [08-Sep-2020 21:17:10] Completed 107 of 160 simulation runs [08-Sep-2020 21:17:15] Completed 108 of 160 simulation runs [08-Sep-2020 21:17:21] Completed 109 of 160 simulation runs [08-Sep-2020 21:17:26] Completed 110 of 160 simulation runs [08-Sep-2020 21:17:32] Completed 111 of 160 simulation runs [08-Sep-2020 21:17:37] Completed 112 of 160 simulation runs [08-Sep-2020 21:17:43] Completed 113 of 160 simulation runs [08-Sep-2020 21:17:48] Completed 114 of 160 simulation runs [08-Sep-2020 21:17:54] Completed 115 of 160 simulation runs [08-Sep-2020 21:18:00] Completed 116 of 160 simulation runs [08-Sep-2020 21:18:05] Completed 117 of 160 simulation runs [08-Sep-2020 21:18:12] Completed 118 of 160 simulation runs [08-Sep-2020 21:18:18] Completed 119 of 160 simulation runs [08-Sep-2020 21:18:24] Completed 120 of 160 simulation runs [08-Sep-2020 21:18:30] Completed 121 of 160 simulation runs [08-Sep-2020 21:18:36] Completed 122 of 160 simulation runs [08-Sep-2020 21:18:42] Completed 123 of 160 simulation runs [08-Sep-2020 21:18:49] Completed 124 of 160 simulation runs [08-Sep-2020 21:18:57] Completed 125 of 160 simulation runs [08-Sep-2020 21:19:03] Completed 126 of 160 simulation runs [08-Sep-2020 21:19:09] Completed 127 of 160 simulation runs [08-Sep-2020 21:19:14] Completed 128 of 160 simulation runs [08-Sep-2020 21:19:20] Completed 129 of 160 simulation runs [08-Sep-2020 21:19:28] Completed 130 of 160 simulation runs [08-Sep-2020 21:19:35] Completed 131 of 160 simulation runs [08-Sep-2020 21:19:41] Completed 132 of 160 simulation runs [08-Sep-2020 21:19:47] Completed 133 of 160 simulation runs [08-Sep-2020 21:19:54] Completed 134 of 160 simulation runs [08-Sep-2020 21:20:03] Completed 135 of 160 simulation runs [08-Sep-2020 21:20:12] Completed 136 of 160 simulation runs [08-Sep-2020 21:20:19] Completed 137 of 160 simulation runs [08-Sep-2020 21:20:27] Completed 138 of 160 simulation runs [08-Sep-2020 21:20:35] Completed 139 of 160 simulation runs [08-Sep-2020 21:20:43] Completed 140 of 160 simulation runs [08-Sep-2020 21:20:50] Completed 141 of 160 simulation runs [08-Sep-2020 21:20:56] Completed 142 of 160 simulation runs [08-Sep-2020 21:21:03] Completed 143 of 160 simulation runs [08-Sep-2020 21:21:13] Completed 144 of 160 simulation runs [08-Sep-2020 21:21:19] Completed 145 of 160 simulation runs [08-Sep-2020 21:21:25] Completed 146 of 160 simulation runs [08-Sep-2020 21:21:32] Completed 147 of 160 simulation runs [08-Sep-2020 21:21:42] Completed 148 of 160 simulation runs [08-Sep-2020 21:21:49] Completed 149 of 160 simulation runs [08-Sep-2020 21:21:57] Completed 150 of 160 simulation runs [08-Sep-2020 21:22:05] Completed 151 of 160 simulation runs [08-Sep-2020 21:22:14] Completed 152 of 160 simulation runs [08-Sep-2020 21:22:22] Completed 153 of 160 simulation runs [08-Sep-2020 21:22:31] Completed 154 of 160 simulation runs [08-Sep-2020 21:22:39] Completed 155 of 160 simulation runs [08-Sep-2020 21:22:47] Completed 156 of 160 simulation runs [08-Sep-2020 21:22:55] Completed 157 of 160 simulation runs [08-Sep-2020 21:23:03] Completed 158 of 160 simulation runs [08-Sep-2020 21:23:12] Completed 159 of 160 simulation runs [08-Sep-2020 21:23:21] Completed 160 of 160 simulation runs [08-Sep-2020 21:23:23] Cleaning up parallel workers...
ens = simulationEnsembleDatastore(fullfile('.','Data'));
Модель сконфигурирована, чтобы регистрировать давление выхода насоса, вывести поток, частоту вращения двигателя и моторный ток.
ens.DataVariables
ans = 8×1 string
"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));
Evaluating tall expression using the Parallel Pool 'local': - Pass 1 of 1: 0% complete Evaluation 0% complete
- Pass 1 of 1: Completed in 11 sec Evaluation completed in 11 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.84932 117.83 18.867 276.49 35.573 7.5235 -0.73064 2.778 13.86 1.1491 35.679 2.2319 42692 1e-09 0.8 0
43.909 0.46288 126.04 18.969 12.417 35.577 7.8766 -0.70939 2.63 13.336 1.1451 35.688 2.3235 42699 1e-09 0.8 0
44.03 66.853 151.07 23.624 230.13 34.357 9.4974 -0.54522 2.8046 15.243 1.1665 34.495 2.4703 41231 1.6e-06 0.71 0
40.699 6.0325 89.768 21.02 252.48 32.617 6.7424 -0.70295 2.7538 12.264 1.1457 32.721 2.1208 39139 4e-07 0.8 6.6667e-05
8.6012 14.366 29.295 9.8881 485.93 18.623 9.325 -0.02553 2.0546 14.315 1.3422 18.871 2.6236 22347 2.8e-06 0.8 0.0006
43.969 33.122 141.94 24.284 154.8 34.682 8.6405 -0.53746 2.7521 15.692 1.1955 34.807 2.3614 41620 1.2e-06 0.8 0
9.5096 25.219 31.388 13.27 65.397 21.352 6.7393 -0.25342 2.4885 13.704 1.2738 21.509 2.1357 25619 2e-06 0.8 0.00046667
26.83 1.114 31.794 15.963 434.6 21.432 3.6724 -0.44013 2.6938 10.304 1.2108 21.517 1.5589 25716 4e-07 0.8 0.00053333
9.0252 16.051 23.779 12.955 302.15 20.006 7.5907 -0.13753 2.2392 14.006 1.2837 20.195 2.3023 24010 2.4e-06 0.8 0.00053333
43.969 12.291 124.8 22.447 478.48 34.986 8.0954 -0.71245 2.9268 15.015 1.1585 35.101 2.2984 41987 8e-07 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.5000
% Compute the accuracy of the classifier at predicting % that there is a fault 1-sum(confdata(2:end,1))/sum(confdata(:))
ans = 0.9375
Исследуйте случаи, где никакой отказ не был предсказан, но отказ действительно существовал. Сначала найдите случаи в данных о валидации, где фактический отказ был блокирующимся отказом, но никакой отказ не был предсказан.
vData = data(cvp.test(1),:); b1 = (actualValue==2) & (predictedValue==0); fData = vData(b1,15:17)
fData=2×3 table
LeakFault BlockingFault BearingFault
_________ _____________ ____________
1e-09 0.74 0
8e-07 0.59 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: 0% complete Evaluation 0% complete
- Pass 1 of 1: Completed in 9.7 sec Evaluation completed in 9.8 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