Используя Simulink, чтобы сгенерировать данные об отказе
Этот пример показывает, как использовать модель Simulink, чтобы сгенерировать отказ и здоровые данные. Отказ и здоровые данные используются, чтобы разработать алгоритм мониторинга состояния. Пример использует систему передачи и моделирует зубной отказ механизма, отказ дрейфа датчика и отказ износа вала.
Системная модель передачи
Модель преобразования регистра передачи использует блоки Simscape Driveline, чтобы смоделировать простую систему передачи. Система передачи состоит из диска крутящего момента, карданного вала, муфты и высоких и низких механизмов, соединенных с выходным валом.
mdl = 'pdmTransmissionCasing';
open_system(mdl)
Система передачи включает датчик вибрации, который контролирует случающиеся колебания. Случающаяся модель переводит вал угловое смещение в линейное смещение на преобразовании регистра. Преобразование регистра моделируется как массовая пружинная система демпфера, и вибрация (заключающий ускорение в корпус) измеряется от преобразования регистра.
open_system([mdl '/Casing'])
Моделирование отказа
Система передачи включает модели отказа для дрейфа датчика вибрации, зубного отказа механизма и износа вала. Отказ дрейфа датчика легко моделируется путем представления смещения в модели датчика. Смещением управляет образцовая переменная SDrift, обратите внимание, что SDrift=0 не подразумевает отказа датчика.
open_system([mdl '/Vibration sensor with drift'])
Отказ износа вала моделируется различной подсистемой. В этом случае варианты подсистемы изменяют затухание вала, но различная подсистема могла использоваться, чтобы полностью изменить реализацию модели вала. Выбранным вариантом управляет модель variable ShaftWear, обратите внимание, что that 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, чтобы запустить симуляции. Функция generateSimulationEnsemble конфигурирует модель, чтобы сохранить записанные данные в файл, использовать формат расписания для журналирования сигнала и хранить объекты Simulink.SimulationInput в сохраненном файле. Функция generateSimulationEnsemble возвращает флаг состояния, указывающий ли симуляция, завершенная успешно.
Код выше создал 110 входных параметров симуляции из значений переменных с координатной сеткой и входные параметры симуляции 98 от значений случайной переменной, дающих 208 общих симуляций. При выполнении этих 208 симуляций в параллельных взятиях, приблизительно десять минут на стандартном рабочем столе и генерируют приблизительно 10 ГБ данных, возможность, чтобы только запустить первые 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...
generateSimulationEnsemble запустил и регистрировал результаты симуляции. Создайте ансамбль симуляции, чтобы обработать и анализировать результаты симуляции с помощью команды simulationEnsembleDatastore.
ens = simulationEnsembleDatastore(fullfile(pwd,'Data'));Обработка результатов симуляции
Команда simulationEnsembleDatastore создала объект ансамбля, который указывает на результаты симуляции. Используйте объект ансамбля подготовить и анализировать данные в каждом члене ансамбля. Списки объектов ансамбля переменные данных в ансамбле и по умолчанию всех переменных выбраны для чтения.
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. 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);
Сигнал 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
Свойство 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Обработанная вибрация и сигналы tacho и метки отказа добавляются к ансамблю, чтобы использоваться позже.
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 управляет, чтобы цикл использовался, чтобы применить 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"];
Поскольку каждый участник в ансамбле вычисляет много времен и основанных на спектре функций. Они включают статистику сигнала, такую как среднее значение сигнала, отклонение, пик к пиковым, нелинейным характеристикам сигнала, таким как аппроксимированная энтропия и экспонента Ляпунова и спектральные функции, такие как пиковая частота времени синхронное среднее значение сигнала вибрации и степень времени синхронный средний сигнал конверта. Функция 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 со всеми функциями, вычисленными выше как предикторы и метка отказа дрейфа датчика (истинное ложное значение) как ответ. Команда 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, с помощью ансамбля симуляции, чтобы очистить данные моделирования и функции извлечения. Извлеченные функции были затем использованы, чтобы создать классификаторы для различных типов отказа.