Использование Simulink для генерации данных о отказе
В этом примере показано, как использовать модель Simulink для генерации отказов и исправных данных. Данные об отказе и исправности используются для разработки алгоритма мониторинга условия. Пример использует систему трансмиссии и моделирует отказ зуба шестерни, отказ дрейфа датчика и отказ износа вала.
Модель трансмиссии
Модель корпуса трансмиссии использует блоки Simscape Driveline, чтобы смоделировать простую систему трансмиссии. Система трансмиссии состоит из привода крутящего момента, приводного вала, сцепления и высоких и низких передач, соединенных с выходом валом.
mdl = 'pdmTransmissionCasing';
open_system(mdl)
Система трансмиссии включает датчик вибрации, который контролирует вибрации корпуса. Модель корпуса перемещает угловое перемещение вала к линейному перемещению на корпусе. Корпус моделируется как система демпфера весовой пружины, и вибрация (ускорение корпуса) измеряется от корпуса.
open_system([mdl '/Casing'])
Моделирование отказа
Система трансмиссии включает модели отказа для дрейфа датчика вибрации, отказа зуба шестерни и износа вала. Отказ дрейфа датчика легко моделируется путем введения смещения в модель датчика. Смещением управляет переменная модели SDrift, обратите внимание, что SDrift=0 не подразумевает отказа датчика.
open_system([mdl '/Vibration sensor with drift'])
Отказ износа вала моделируется альтернативной подсистемой. В этом случае варианты подсистемы изменяют демпфирование вала, но альтернативная подсистема может использоваться, чтобы полностью изменить реализацию модели вала. Выбранный вариант управляется переменной модели ShaftWear, обратите внимание, что ShaftWear=0 означает отсутствие отказа вала.
open_system([mdl '/Shaft'])
open_system([mdl,'/Shaft/Output Shaft'])
Отказ зубчатого колеса моделируется инжекцией нарушения порядка при фиксированном положении во вращении приводного вала. Положение вала измеряется в радианах, и когда положение вала находится внутри небольшого окна около 0, к валу прикладывается нарушение порядка. Величина нарушения порядка определяется переменной модели ToothFaultGain, обратите внимание, что ToothFaultGain=0 не подразумевает отказа зубчатого колеса.
Симуляция отказа и исправность данных
Модель трансмиссии сконфигурирована с переменными, которые контролируют присутствие и серьезность трех различных типов отказов, дрейфа датчика, зуба передачи и износа вала. Варьируя переменные модели, SDrift, ToothFaultGain, и ShaftWearможно создать данные о вибрации для различных типов отказов. Используйте массив Simulink.SimulationInput объекты для определения ряда различных сценариев симуляции. Например, выберите массив значений для каждой из переменных модели, а затем используйте ndgrid функция для создания Simulink.SimulationInput для каждой комбинации значений переменных модели.
toothFaultArray = -2:2/10:0; % Tooth fault gain values sensorDriftArray = -1:0.5:1; % Sensor drift offset values shaftWearArray = [0 -1]; % Variants available for drive shaft conditions % Create an n-dimensional array with combinations of all values [toothFaultValues,sensorDriftValues,shaftWearValues] = ... ndgrid(toothFaultArray,sensorDriftArray,shaftWearArray); for ct = numel(toothFaultValues):-1:1 % Create a Simulink.SimulationInput for each combination of values siminput = Simulink.SimulationInput(mdl); % Modify model parameters siminput = setVariable(siminput,'ToothFaultGain',toothFaultValues(ct)); siminput = setVariable(siminput,'SDrift',sensorDriftValues(ct)); siminput = setVariable(siminput,'ShaftWear',shaftWearValues(ct)); % Collect the simulation input in an array gridSimulationInput(ct) = siminput; end
Точно так же создайте комбинации случайных значений для каждой переменной модели. Обязательно включите 0 значение так, что существуют комбинации, где представлены только подмножества из трех типов отказов.
rng('default'); % Reset random seed for reproducibility toothFaultArray = [0 -rand(1,6)]; % Tooth fault gain values sensorDriftArray = [0 randn(1,6)/8]; % Sensor drift offset values shaftWearArray = [0 -1]; % Variants available for drive shaft conditions %Create an n-dimensional array with combinations of all values [toothFaultValues,sensorDriftValues,shaftWearValues] = ... ndgrid(toothFaultArray,sensorDriftArray,shaftWearArray); for ct=numel(toothFaultValues):-1:1 % Create a Simulink.SimulationInput for each combination of values siminput = Simulink.SimulationInput(mdl); % Modify model parameters siminput = setVariable(siminput,'ToothFaultGain',toothFaultValues(ct)); siminput = setVariable(siminput,'SDrift',sensorDriftValues(ct)); siminput = setVariable(siminput,'ShaftWear',shaftWearValues(ct)); % Collect the simulation input in an array randomSimulationInput(ct) = siminput; end
С массивом Simulink.SimulationInput объекты, заданные как, используют generateSimulationEnsemble функция для запуска симуляций. The generateSimulationEnsemble функция конфигурирует модель, чтобы сохранить записанные данные в файл, использовать формат timetable для логгирования сигналов и хранить Simulink.SimulationInput объекты в сохраненном файле. The generateSimulationEnsemble функция возвращает флаг состояния, указывающий, завершена ли симуляция успешно.
Приведенный выше код создал 110 симуляции входов из сетчатых значений переменных и 98 входы симуляции из значений случайных переменных, дающих 208 суммарных симуляций. Выполнение этих 208 параллельных симуляций может занять аж два часа или более на стандартном рабочем столе и генерирует около 10GB данных. Опция запуска только первых 10 симуляций предусмотрена для удобства.
% Run the simulations and create an ensemble to manage the simulation results if ~exist(fullfile(pwd,'Data'),'dir') mkdir(fullfile(pwd,'Data')) % Create directory to store results end runAll = true; if runAll [ok,e] = generateSimulationEnsemble([gridSimulationInput, randomSimulationInput], ... fullfile(pwd,'Data'),'UseParallel', true); else [ok,e] = generateSimulationEnsemble(gridSimulationInput(1:10), fullfile(pwd,'Data')); %#ok<*UNRCH> end
[28-Feb-2018 13:17:31] Checking for availability of parallel pool... [28-Feb-2018 13:17:37] Loading Simulink on parallel workers... Analyzing and transferring files to the workers ...done. [28-Feb-2018 13:17:56] Configuring simulation cache folder on parallel workers... [28-Feb-2018 13:17:57] Running SetupFcn on parallel workers... [28-Feb-2018 13:18:01] Loading model on parallel workers... [28-Feb-2018 13:18:02] Transferring base workspace variables used in the model to parallel workers... [28-Feb-2018 13:18:07] Running simulations... [28-Feb-2018 13:18:29] Completed 1 of 208 simulation runs [28-Feb-2018 13:18:29] Completed 2 of 208 simulation runs [28-Feb-2018 13:18:30] Completed 3 of 208 simulation runs ... [28-Feb-2018 13:24:22] Completed 204 of 208 simulation runs [28-Feb-2018 13:24:28] Completed 205 of 208 simulation runs [28-Feb-2018 13:24:28] Completed 206 of 208 simulation runs [28-Feb-2018 13:24:28] Completed 207 of 208 simulation runs [28-Feb-2018 13:24:29] Completed 208 of 208 simulation runs [28-Feb-2018 13:24:33] Cleaning up parallel workers...
The generateSimulationEnsemble запустил и зарегистрировал результаты симуляции. Создайте симуляционный ансамбль для обработки и анализа результатов симуляции с помощью simulationEnsembleDatastore команда.
ens = simulationEnsembleDatastore(fullfile(pwd,'Data'));Обработка результатов симуляции
The simulationEnsembleDatastore команда создала объект ансамбля, который указывает на результаты симуляции. Используйте объект ансамбля для подготовки и анализа данных в каждом представителе ансамбля. Объект ensemble перечисляет переменные данных в ансамбле, и по умолчанию для чтения выбираются все переменные.
ens
ens =
simulationEnsembleDatastore with properties:
DataVariables: [6×1 string]
IndependentVariables: [0×0 string]
ConditionVariables: [0×0 string]
SelectedVariables: [6×1 string]
NumMembers: 208
LastMemberRead: [0×0 string]
ens.SelectedVariables
ans = 6×1 string array
"SimulationInput"
"SimulationMetadata"
"Tacho"
"Vibration"
"xFinal"
"xout"
Для анализа только читайте Vibration и Tacho сигналы и Simulink.SimulationInput. The Simulink.SimulationInput имеет значения переменных модели, используемые для симуляции, и используется для создания меток отказа для представителей ансамбля. Используйте ансамбль read команда для получения данных о представителе ансамбля.
ens.SelectedVariables = ["Vibration" "Tacho" "SimulationInput"]; data = read(ens)
data=1×3 table
Vibration Tacho SimulationInput
___________________ ___________________ ______________________________
[40272×1 timetable] [40272×1 timetable] [1×1 Simulink.SimulationInput]
Извлеките сигнал вибрации из возвращенных данных и постройте график.
vibration = data.Vibration{1};
plot(vibration.Time,vibration.Data)
title('Vibration')
ylabel('Acceleration')
Первые 10 секунд симуляции содержат данные, где система передачи запускается; для анализа отменить эти данные.
idx = vibration.Time >= seconds(10); vibration = vibration(idx,:); vibration.Time = vibration.Time - vibration.Time(1);
The Tacho сигнал содержит импульсы для вращений привода и валов нагрузки, но анализ, а в частности синхронное среднее по времени, требует времени вращений вала. Следующий код отбрасывает первые 10 секунд Tacho данные и находит время вращения вала в tachoPulses.
tacho = data.Tacho{1};
idx = tacho.Time >= seconds(10);
tacho = tacho(idx,:);
plot(tacho.Time,tacho.Data)
title('Tacho pulses')
legend('Drive shaft','Load shaft') % Load shaft rotates more slowly than drive shaft
idx = diff(tacho.Data(:,2)) > 0.5; tachoPulses = tacho.Time(find(idx)+1)-tacho.Time(1)
tachoPulses = 8×1 duration array
2.8543 sec
6.6508 sec
10.447 sec
14.244 sec
18.04 sec
21.837 sec
25.634 sec
29.43 sec
The Simulink.SimulationInput.Variables свойство содержит значения параметров отказа, используемых для симуляции, эти значения позволяют нам создавать метки отказа для каждого представителя ансамбля.
vars = data.SimulationInput{1}.Variables;
idx = strcmp({vars.Name},'SDrift');
if any(idx)
sF = abs(vars(idx).Value) > 0.01; % Small drift values are not faults
else
sF = false;
end
idx = strcmp({vars.Name},'ShaftWear');
if any(idx)
sV = vars(idx).Value < 0;
else
sV = false;
end
if any(idx)
idx = strcmp({vars.Name},'ToothFaultGain');
sT = abs(vars(idx).Value) < 0.1; % Small tooth fault values are not faults
else
sT = false
end
faultCode = sF + 2*sV + 4*sT; % A fault code to capture different fault conditionsОбработанные сигналы вибрации и тахометра и метки отказа добавляются в ансамбль, который будет использоваться позже.
sdata = table({vibration},{tachoPulses},sF,sV,sT,faultCode, ...
'VariableNames',{'Vibration','TachoPulses','SensorDrift','ShaftWear','ToothFault','FaultCode'}) sdata=1×6 table
Vibration TachoPulses SensorDrift ShaftWear ToothFault FaultCode
___________________ ______________ ___________ _________ __________ _________
[30106×1 timetable] [8×1 duration] true false false 1
ens.DataVariables = [ens.DataVariables; "TachoPulses"];Ансамбль ConditionVariables свойство может использоваться, чтобы идентифицировать переменные в ансамбле, которые содержат данные о условии или метке отказа. Установите свойство таким образом, чтобы оно содержало вновь созданные метки отказа.
ens.ConditionVariables = ["SensorDrift","ShaftWear","ToothFault","FaultCode"];
Приведенный выше код использовался для обработки одного представителя ансамбля. Для обработки всех представителей ансамбля приведенный выше код был преобразован в функцию prepareData и использование ансамбля hasdata команда a loop используется для применения prepareData всем представителям ансамбля. Представители ансамбля могут быть обработаны параллельно путем разбиения ансамбля и параллельной обработки ансамблевых перегородок.
reset(ens) runLocal = false; if runLocal % Process each member in the ensemble while hasdata(ens) data = read(ens); addData = prepareData(data); writeToLastMemberRead(ens,addData) end else % Split the ensemble into partitions and process each partition in parallel n = numpartitions(ens,gcp); parfor ct = 1:n subens = partition(ens,n,ct); while hasdata(subens) data = read(subens); addData = prepareData(data); writeToLastMemberRead(subens,addData) end end end
Постройте график сигнала вибрации каждого 10-ого представителя ансамбля с помощью hasdata и read команды для извлечения сигнала вибрации.
reset(ens) ens.SelectedVariables = "Vibration"; figure, ct = 1; while hasdata(ens) data = read(ens); if mod(ct,10) == 0 vibration = data.Vibration{1}; plot(vibration.Time,vibration.Data) hold on end ct = ct + 1; end hold off title('Vibration signals') ylabel('Acceleration')
Анализ данных моделирования
Теперь, когда данные были очищены и предварительно обработаны, данные готовы к извлечению признаков, чтобы определить функции, используемые для классификации различных типов отказов. Сначала сконфигурируйте ансамбль так, чтобы он возвращал только обработанные данные.
ens.SelectedVariables = ["Vibration","TachoPulses"];
Для каждого представителя в ансамбле вычислите количество временных и спектральных функций. Они включают статистику сигналов, такую как среднее значение сигнала, отклонение, пик-пик, нелинейные характеристики сигнала, такие как приблизительная энтропия и экспонента Ляпунова, и спектральные функции, такие как пиковая частота синхронного временного среднего сигнала вибрации и степени синхронного среднего огибающего сигнала. The analyzeData функция содержит полный список извлечённых функций. В качестве примера рассмотрим вычисление спектра синхронного усредненного сигнала вибрации во времени.
reset(ens) data = read(ens)
data=1×2 table
Vibration TachoPulses
___________________ ______________
[30106×1 timetable] [8×1 duration]
vibration = data.Vibration{1};
% Interpolate the Vibration signal onto periodic time base suitable for fft analysis
np = 2^floor(log(height(vibration))/log(2));
dt = vibration.Time(end)/(np-1);
tv = 0:dt:vibration.Time(end);
y = retime(vibration,tv,'linear');
% Compute the time synchronous average of the vibration signal
tp = seconds(data.TachoPulses{1});
vibrationTSA = tsa(y,tp);
figure
plot(vibrationTSA.ttTime,vibrationTSA.tsa)
title('Vibration time synchronous average')
ylabel('Acceleration')
% Compute the spectrum of the time synchronous average np = numel(vibrationTSA); f = fft(vibrationTSA.tsa.*hamming(np))/np; frTSA = f(1:floor(np/2)+1); % TSA frequency response wTSA = (0:np/2)/np*(2*pi/seconds(dt)); % TSA spectrum frequencies mTSA = abs(frTSA); % TSA spectrum magnitudes figure semilogx(wTSA,20*log10(mTSA)) title('Vibration spectrum') xlabel('rad/s')
Частота, которая соответствует пиковой величине, может оказаться функцией, которая полезна для классификации различных типов отказов. Приведенный ниже код вычисляет функции, упомянутые выше, для всех представителей ансамбля (выполнение этого анализа может занять до часа на стандартном рабочем столе. Необязательный код для параллельного запуска анализа с помощью ансамбля partition команда предусмотрена.) Имена элементов добавляются к свойству переменных данных ансамбля перед вызовом writeToLastMemberRead, чтобы добавить вычисленные функции к каждому представителю ансамбля.
reset(ens) ens.DataVariables = [ens.DataVariables; ... "SigMean";"SigMedian";"SigRMS";"SigVar";"SigPeak";"SigPeak2Peak";"SigSkewness"; ... "SigKurtosis";"SigCrestFactor";"SigMAD";"SigRangeCumSum";"SigCorrDimension";"SigApproxEntropy"; ... "SigLyapExponent";"PeakFreq";"HighFreqPower";"EnvPower";"PeakSpecKurtosis"]; if runLocal while hasdata(ens) data = read(ens); addData = analyzeData(data); writeToLastMemberRead(ens,addData); end else % Split the ensemble into partitions and analyze each partition in parallel n = numpartitions(ens,gcp); parfor ct = 1:n subens = partition(ens,n,ct); while hasdata(subens) data = read(subens); addData = analyzeData(data); writeToLastMemberRead(subens,addData) end end end
Выбор функций для классификации отказов
Функции, вычисленные выше, используются для создания классификатора для классификации различных условий отказа. Сначала сконфигурируйте ансамбль, чтобы считать только производные функции и метки отказа.
featureVariables = analyzeData('GetFeatureNames');
ens.DataVariables = [ens.DataVariables; featureVariables];
ens.SelectedVariables = [featureVariables; ens.ConditionVariables];
reset(ens)Соберите данные функций для всех представителей ансамбля в одну таблицу.
featureData = gather(tall(ens))
Evaluating tall expression using the Parallel Pool 'local': - Pass 1 of 1: Completed in 1.8 min Evaluation completed in 1.8167 min
featureData=208×22 table
SigMean SigMedian SigRMS SigVar SigPeak SigPeak2Peak SigSkewness SigKurtosis SigCrestFactor SigMAD SigRangeCumSum SigCorrDimension SigApproxEntropy SigLyapExponent PeakFreq HighFreqPower EnvPower PeakSpecKurtosis SensorDrift ShaftWear ToothFault FaultCode
________ _________ _______ _______ _______ ____________ ___________ ___________ ______________ _______ ______________ ________________ ________________ _______________ ________ _____________ __________ ________________ ___________ _________ __________ _________
-0.94876 -0.9722 1.3726 0.98387 0.81571 3.6314 -0.041525 2.2666 2.0514 0.8081 28562 1.1427 0.031601 79.531 0 6.7528e-06 3.2349e-07 162.13 true false false 1
-0.97537 -0.98958 1.3937 0.99105 0.81571 3.6314 -0.023777 2.2598 2.0203 0.81017 29418 1.1362 0.03786 70.339 0 5.0788e-08 9.1568e-08 226.12 true false false 1
1.0502 1.0267 1.4449 0.98491 2.8157 3.6314 -0.04162 2.2658 1.9487 0.80853 31710 1.1479 0.031586 125.18 0 6.7416e-06 3.1343e-07 162.13 true true false 3
1.0227 1.0045 1.4288 0.99553 2.8157 3.6314 -0.016356 2.2483 1.9707 0.81324 30984 1.1472 0.032109 112.52 0 4.9934e-06 2.5787e-07 162.13 true true false 3
1.0123 1.0024 1.4202 0.99233 2.8157 3.6314 -0.014701 2.2542 1.9826 0.81156 30661 1.147 0.032891 109.02 0 3.619e-06 2.2397e-07 230.39 true true false 3
1.0275 1.0102 1.4338 1.0001 2.8157 3.6314 -0.02659 2.2439 1.9638 0.81589 31102 1.0975 0.033449 64.499 0 2.5493e-06 1.9224e-07 230.39 true true false 3
1.0464 1.0275 1.4477 1.0011 2.8157 3.6314 -0.042849 2.2455 1.9449 0.81595 31665 1.1417 0.034182 98.759 0 1.7313e-06 1.6263e-07 230.39 true true false 3
1.0459 1.0257 1.4402 0.98047 2.8157 3.6314 -0.035405 2.2757 1.955 0.80583 31554 1.1345 0.035323 44.304 0 1.1115e-06 1.2807e-07 230.39 true true false 3
1.0263 1.0068 1.4271 0.98341 2.8157 3.6314 -0.0165 2.2726 1.973 0.80624 30951 1.1515 0.03592 125.28 0 6.5947e-07 1.208e-07 162.13 true true false 3
1.0103 1.0014 1.4183 0.99091 2.8157 3.6314 -0.011667 2.2614 1.9853 0.80987 30465 1.0619 0.036514 17.093 0 5.2297e-07 1.0704e-07 230.39 true true false 3
1.0129 1.0023 1.419 0.98764 2.8157 3.6314 -0.010969 2.266 1.9843 0.80866 30523 1.1371 0.037234 84.568 0 2.1605e-07 8.774e-08 230.39 true true false 3
1.0251 1.0104 1.4291 0.9917 2.8157 3.6314 -0.023609 2.2588 1.9702 0.81048 30896 1.1261 0.037836 98.463 0 5.0275e-08 9.0495e-08 230.39 true true false 3
-0.97301 -0.99243 1.3928 0.99326 0.81571 3.6314 -0.012569 2.2589 2.0216 0.81095 29351 1.1277 0.038507 42.887 0 1.1383e-11 8.3005e-08 230.39 true false true 5
1.0277 1.0084 1.4315 0.99314 2.8157 3.6314 -0.013542 2.2598 1.9669 0.81084 30963 1.1486 0.038524 99.426 0 1.1346e-11 8.353e-08 226.12 true true true 7
0.026994 0.0075709 0.99697 0.99326 1.8157 3.6314 -0.012569 2.2589 1.8212 0.81095 1083.8 1.1277 0.038507 44.448 9.998 4.9172e-12 8.3005e-08 230.39 false false true 4
0.026943 0.0084639 0.99297 0.98531 1.8157 3.6314 -0.018182 2.2686 1.8286 0.80732 1466.6 1.1368 0.035822 93.95 1.8618 6.8872e-07 1.2185e-07 230.39 false false false 0
⋮
Учитывайте отказ дрейфа датчика. Используйте fscnca команда со всеми функциями, вычисленными выше как предикторы, и метка отказа дрейфа датчика (истинное ложное значение) как ответ. The fscnca команда возвращает веса для каждой функции и признаки с более высокими весами имеют более большое значение в прогнозировании отклика. Для отказа дрейфа датчика веса указывают, что две функции являются значительными предикторами (совокупная сумма сигнала области значений и пиковая частота спектрального куртоза), а остальные функции оказывают незначительное влияние.
idxResponse = strcmp(featureData.Properties.VariableNames,'SensorDrift'); idxLastFeature = find(idxResponse)-1; % Index of last feature to use as a potential predictor featureAnalysis = fscnca(featureData{:,1:idxLastFeature},featureData.SensorDrift); featureAnalysis.FeatureWeights
ans = 18×1
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
⋮
idxSelectedFeature = featureAnalysis.FeatureWeights > 0.1; classifySD = [featureData(:,idxSelectedFeature), featureData(:,idxResponse)]
classifySD=208×3 table
SigRangeCumSum PeakSpecKurtosis SensorDrift
______________ ________________ ___________
28562 162.13 true
29418 226.12 true
31710 162.13 true
30984 162.13 true
30661 230.39 true
31102 230.39 true
31665 230.39 true
31554 230.39 true
30951 162.13 true
30465 230.39 true
30523 230.39 true
30896 230.39 true
29351 230.39 true
30963 226.12 true
1083.8 230.39 false
1466.6 230.39 false
⋮
Сгруппированная гистограмма совокупной суммы области значений дает нам представление о том, почему эта функция является значимым предиктором отказа дрейфа датчика.
figure histogram(classifySD.SigRangeCumSum(classifySD.SensorDrift),'BinWidth',5e3) xlabel('Signal cumulative sum range') ylabel('Count') hold on histogram(classifySD.SigRangeCumSum(~classifySD.SensorDrift),'BinWidth',5e3) hold off legend('Sensor drift fault','No sensor drift fault')
График гистограммы показывает, что совокупный суммарный диапазон сигнала является хорошим показателем для обнаружения сбоев дрейфа датчика, хотя, вероятно, необходима дополнительная функция, поскольку могут быть ложные срабатывания, когда совокупный суммарный диапазон сигнала ниже 0,5, если только совокупный суммарный диапазон сигнала используется для классификации дрейфа датчика.
Учитывайте отказ износа вала. В этом случае fscnca функция указывает, что существует 3 признака, которые являются значимыми предикторами отказа (сигнал Ляпунова экспонента, пиковая частота и пиковая частота в спектральном куртозе), выберите их, чтобы классифицировать отказ износа вала.
idxResponse = strcmp(featureData.Properties.VariableNames,'ShaftWear');
featureAnalysis = fscnca(featureData{:,1:idxLastFeature},featureData.ShaftWear);
featureAnalysis.FeatureWeightsans = 18×1
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
⋮
idxSelectedFeature = featureAnalysis.FeatureWeights > 0.1; classifySW = [featureData(:,idxSelectedFeature), featureData(:,idxResponse)]
classifySW=208×4 table
SigLyapExponent PeakFreq PeakSpecKurtosis ShaftWear
_______________ ________ ________________ _________
79.531 0 162.13 false
70.339 0 226.12 false
125.18 0 162.13 true
112.52 0 162.13 true
109.02 0 230.39 true
64.499 0 230.39 true
98.759 0 230.39 true
44.304 0 230.39 true
125.28 0 162.13 true
17.093 0 230.39 true
84.568 0 230.39 true
98.463 0 230.39 true
42.887 0 230.39 false
99.426 0 226.12 true
44.448 9.998 230.39 false
93.95 1.8618 230.39 false
⋮
Сгруппированная гистограмма для экспоненты сигнала Ляпунова показывает, почему одна эта функция не является хорошим предиктором.
figure histogram(classifySW.SigLyapExponent(classifySW.ShaftWear)) xlabel('Signal lyapunov exponent') ylabel('Count') hold on histogram(classifySW.SigLyapExponent(~classifySW.ShaftWear)) hold off legend('Shaft wear fault','No shaft wear fault')
Выбор признаков износа вала указывает, что для классификации отказа износа вала необходимо несколько функций, сгруппированная гистограмма подтверждает это как наиболее значимейшая функция (в этом случае показатель Ляпунова) имеет аналогичное распределение как для сценариев отказа, так и для сценариев без отказа, указывающих, что для правильной классификации этого отказа необходимо больше функций.
Наконец, рассмотрим отказ зуба, fscnca функция указывает, что существует 3 первичных признака, которые являются значительными предикторами отказа (совокупная сумма сигнала области значений, экспонента Ляпунова сигнала и пиковая частота в спектральном куртозе). Выбор этих 3 функций для классификации отказа зуба приводит к классификатору, который имеет низкую эффективность. Вместо этого используйте 6 наиболее важных функций.
idxResponse = strcmp(featureData.Properties.VariableNames,'ToothFault');
featureAnalysis = fscnca(featureData{:,1:idxLastFeature},featureData.ToothFault);
[~,idxSelectedFeature] = sort(featureAnalysis.FeatureWeights);
classifyTF = [featureData(:,idxSelectedFeature(end-5:end)), featureData(:,idxResponse)]classifyTF=208×7 table
PeakFreq SigPeak SigCorrDimension SigLyapExponent PeakSpecKurtosis SigRangeCumSum ToothFault
________ _______ ________________ _______________ ________________ ______________ __________
0 0.81571 1.1427 79.531 162.13 28562 false
0 0.81571 1.1362 70.339 226.12 29418 false
0 2.8157 1.1479 125.18 162.13 31710 false
0 2.8157 1.1472 112.52 162.13 30984 false
0 2.8157 1.147 109.02 230.39 30661 false
0 2.8157 1.0975 64.499 230.39 31102 false
0 2.8157 1.1417 98.759 230.39 31665 false
0 2.8157 1.1345 44.304 230.39 31554 false
0 2.8157 1.1515 125.28 162.13 30951 false
0 2.8157 1.0619 17.093 230.39 30465 false
0 2.8157 1.1371 84.568 230.39 30523 false
0 2.8157 1.1261 98.463 230.39 30896 false
0 0.81571 1.1277 42.887 230.39 29351 true
0 2.8157 1.1486 99.426 226.12 30963 true
9.998 1.8157 1.1277 44.448 230.39 1083.8 true
1.8618 1.8157 1.1368 93.95 230.39 1466.6 false
⋮
figure histogram(classifyTF.SigRangeCumSum(classifyTF.ToothFault)) xlabel('Signal cumulative sum range') ylabel('Count') hold on histogram(classifyTF.SigRangeCumSum(~classifyTF.ToothFault)) hold off legend('Gear tooth fault','No gear tooth fault')
Использование вышеописанных результатов полинома svm для классификации отказов зуба шестерни. Разделите таблицу функций на представители, которые используются для обучения, и представители для проверки и валидации. Используйте представители, чтобы создать классификатор svm с помощью fitcsvm команда.
rng('default') % For reproducibility cvp = cvpartition(size(classifyTF,1),'KFold',5); % Randomly partition the data for training and validation classifierTF = fitcsvm(... classifyTF(cvp.training(1),1:end-1), ... classifyTF(cvp.training(1),end), ... 'KernelFunction','polynomial', ... 'PolynomialOrder',2, ... 'KernelScale','auto', ... 'BoxConstraint',1, ... 'Standardize',true, ... 'ClassNames',[false; true]);
Используйте классификатор, чтобы классифицировать тестовые точки с помощью predict и проверяйте эффективность предсказаний, используя матрицу неточностей.
% Use the classifier on the test validation data to evaluate performance actualValue = classifyTF{cvp.test(1),end}; predictedValue = predict(classifierTF, classifyTF(cvp.test(1),1:end-1)); % Check performance by computing and plotting the confusion matrix confdata = confusionmat(actualValue,predictedValue); h = heatmap(confdata, ... 'YLabel', 'Actual gear tooth fault', ... 'YDisplayLabels', {'False','True'}, ... 'XLabel', 'Predicted gear tooth fault', ... 'XDisplayLabels', {'False','True'}, ... 'ColorbarVisible','off');
Матрица неточностей указывает, что классификатор правильно классифицирует все условия, не связанные с отказом, но неправильно классифицирует одно ожидаемое состояние отказа как не являющееся отказом. Увеличение количества функций, используемых в классификаторе, может способствовать дальнейшему повышению эффективности.
Сводные данные
Этот пример прошел по рабочему процессу для генерации данных о отказе из Simulink, используя ансамбль симуляции, чтобы очистить данные моделирования и извлечь функции. Извлеченные функции затем использовались для создания классификаторов для различных типов отказов.