Simulink를 사용하여 결함 데이터 생성하기

다음 제품이 필요합니다.

이 예제에서는 Simulink® 모델을 사용하여 결함 데이터와 정상 데이터를 생성하는 방법을 보여줍니다. 결함 데이터와 정상 데이터는 상태 모니터링 알고리즘을 개발하는 데 사용됩니다. 이 예제에서는 변속기 시스템을 사용하여 톱니 결함, 센서 드리프트 결함, 축 마모 결함이 있는 기어를 모델링합니다.

변속기 시스템 모델

변속기 케이싱 모델은 Simscape™ Driveline™ 블록을 사용하여 간단한 변속기 시스템을 모델링합니다. 변속기 시스템은 토크 구동, 구동 축, 클러치, 그리고 출력 축에 연결된 고속 및 저속 기어로 구성됩니다.

mdl = 'pdmTransmissionCasing';
open_system(mdl)

변속기 시스템에는 케이싱 진동을 모니터링하는 진동 센서가 있습니다. 케이싱 모델은 축의 각변위를 케이싱의 선형 변위로 변환합니다. 케이싱은 질량-스프링-댐퍼 시스템으로 모델링되고, 진동(케이싱 가속도)은 케이싱에서 측정됩니다.

open_system([mdl '/Casing'])

결함 모델링

변속기 시스템은 진동 센서 드리프트, 톱니 결함 및 축 마모에 대한 결함 모델을 포함합니다. 센서 드리프트는 센서 모델에 오프셋을 추가하여 쉽게 모델링됩니다. 오프셋은 모델 변수 SDrift에 의해 제어됩니다. 여기서 SDrift=0은 센서 결함이 없음을 나타냅니다.

open_system([mdl '/Vibration sensor with drift'])

축 마모 결함은 Variant Subsystem으로 모델링됩니다. 이 사례에서는 Variant Subsystem이 축 감쇠만 변경하지만, Variant Subsystem을 사용해 축 모델 구현 자체를 완전히 변경할 수도 있습니다. 선택된 Variant는 모델 변수 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개의 시뮬레이션을 병렬로 실행하면 표준 데스크탑에서 2시간 이상 소요될 수 있으며 약 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, 'ShowProgress', false);
else
    [ok,e] = generateSimulationEnsemble(gridSimulationInput(1:10), fullfile(pwd,'Data'), 'ShowProgress', false); %#ok<*UNRCH>
end

Starting parallel pool (parpool) using the 'Processes' profile ...
Preserving jobs with IDs: 1 because they contain crash dump files.
You can use 'delete(myCluster.Jobs)' to remove all jobs created with profile Processes. To create 'myCluster' use 'myCluster = parcluster('Processes')'.
Connected to parallel pool with 6 workers.

generateSimulationEnsemble이 실행되어 시뮬레이션 결과를 기록했습니다. simulationEnsembleDatastore 명령을 사용하여 시뮬레이션 결과를 처리하고 분석하기 위한 시뮬레이션 앙상블을 만듭니다.

ens = simulationEnsembleDatastore(fullfile(pwd,'Data'));

시뮬레이션 결과의 처리

simulationEnsembleDatastore 명령은 시뮬레이션 결과를 가리키는 ensemble 객체를 만들었습니다. 이 ensemble 객체를 사용하여 앙상블의 각 멤버에서 데이터를 준비하고 분석합니다. ensemble 객체는 앙상블에 있는 데이터 변수를 나열하며, 기본적으로 모든 변수가 읽어 올 대상으로 선택됩니다.

ens

ens = 
  simulationEnsembleDatastore with properties:

           DataVariables: [6×1 string]
    IndependentVariables: [0×0 string]
      ConditionVariables: [0×0 string]
       SelectedVariables: [6×1 string]
                ReadSize: 1
              NumMembers: 208
          LastMemberRead: [0×0 string]
                   Files: [208×1 string]

ens.SelectedVariables

ans = 6×1 string
    "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 신호는 구동 축과 부하 축의 회전에 대한 펄스를 포함하지만 분석에서는, 특히 시간 동기 평균화에서는 축 회전의 시간이 필요합니다. 다음 코드는 Tacho 데이터의 처음 10초를 버리고 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
   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

hasdata 명령과 read 명령을 사용하여 진동 신호를 추출하고 앙상블의 10번째 멤버마다 진동 신호를 플로팅합니다.

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.Time,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')

피크 크기에 대응하는 주파수는 여러 결함 유형을 분류하는 데 유용한 특징이 될 수 있습니다. 아래 코드는 모든 앙상블 멤버에 대해 위에서 언급한 특징을 계산합니다. (이 분석을 실행하면 표준 데스크탑에서 최대 1시간이 소요될 수 있습니다. 앙상블 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 'Processes':
- Pass 1 of 1: Completed in 14 sec
Evaluation completed in 14 sec

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.031581            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.037835            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.031565            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.032088            112.52               0      4.9939e-06      2.6719e-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.03287            109.02               0      3.6182e-06      2.1964e-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.033427            64.498               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.034159            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.0353            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.035897            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.036489            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.037209             84.57               0      2.1605e-07      8.5562e-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.037811            98.462               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.038481            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.038499            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.038481            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.035799             93.95          1.8628      6.8862e-07      1.1702e-07         230.39            false         false        false           0    
      ⋮

센서 드리프트 결함을 살펴봅니다. 위에서 계산한 모든 특징을 예측 변수로 지정하고 센서 드리프트 결함 레이블(true/false 값)을 응답 변수로 지정하여 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

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

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.498               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.57               0          230.39           true   
        98.462               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.8628          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
    SigPeak    PeakFreq    SigCorrDimension    SigLyapExponent    PeakSpecKurtosis    SigRangeCumSum    ToothFault
    _______    ________    ________________    _______________    ________________    ______________    __________

    0.81571          0          1.1427             79.531              162.13              28562          false   
    0.81571          0          1.1362             70.339              226.12              29418          false   
     2.8157          0          1.1479             125.18              162.13              31710          false   
     2.8157          0          1.1472             112.52              162.13              30984          false   
     2.8157          0           1.147             109.02              230.39              30661          false   
     2.8157          0          1.0975             64.498              230.39              31102          false   
     2.8157          0          1.1417             98.759              230.39              31665          false   
     2.8157          0          1.1345             44.304              230.39              31554          false   
     2.8157          0          1.1515             125.28              162.13              30951          false   
     2.8157          0          1.0619             17.093              230.39              30465          false   
     2.8157          0          1.1371              84.57              230.39              30523          false   
     2.8157          0          1.1261             98.462              230.39              30896          false   
    0.81571          0          1.1277             42.887              230.39              29351          true    
     2.8157          0          1.1486             99.426              226.12              30963          true    
     1.8157      9.998          1.1277             44.448              230.39             1083.8          true    
     1.8157     1.8628          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이 생성됩니다. 특징 테이블을 훈련용으로 사용되는 멤버와 테스트 및 검증용으로 사용되는 멤버로 분할합니다. fitcsvm 명령을 사용하여 훈련 멤버로 svm 분류기를 만듭니다.

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');