В этом примере показано, как использовать модель Simulink для создания данных об отказах и работоспособности. Данные о неисправности и исправности используются для разработки алгоритма мониторинга состояния. В примере используется трансмиссионная система и моделируется неисправность зубьев шестерен, неисправность дрейфа датчика и неисправность износа вала.
Модель корпуса трансмиссии использует блоки привода Simscape для моделирования простой системы трансмиссии. Трансмиссионная система состоит из привода крутящего момента, ведущего вала, муфты, а также высокой и низкой шестерен, соединенных с выходным валом.
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 для выполнения моделирования. generateSimulationEnsemble функция конфигурирует модель для сохранения зарегистрированных данных в файл, использования формата расписания для регистрации сигнала и сохранения Simulink.SimulationInput объектов в сохраненном файле. 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...
generateSimulationEnsemble выполнил и зарегистрировал результаты моделирования. Создайте ансамбль моделирования для обработки и анализа результатов моделирования с помощью simulationEnsembleDatastore команда.
ens = simulationEnsembleDatastore(fullfile(pwd,'Data')); 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. 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Обработанные сигналы вибрации и тахо-сигналы, а также метки неисправностей добавляются в ансамбль для последующего использования.
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 с использованием ансамбля моделирования для очистки данных моделирования и извлечения элементов. Извлеченные элементы затем использовались для построения классификаторов для различных типов отказов.