Используя 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.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;
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, с помощью ансамбля симуляции, чтобы очистить данные моделирования и функции извлечения. Извлеченные функции были затем использованы, чтобы создать классификаторы для различных типов отказа.

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

Похожие темы

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