ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

Simulink を使用した故障データの生成

この例では、Simulink モデルを使用して健全な状態のデータと故障状態のデータを生成する方法を説明します。故障状態のデータと健全状態のデータを使用して、条件監視のアルゴリズムを開発します。例ではトランスミッション システムを使用してギア歯、センサー ドリフト、およびシャフト摩耗の故障をモデル化します。

トランスミッション システム モデル

トランスミッション ケーシング モデルでは、Simscape Driveline ブロックを使用して単純なトランスミッション システムをモデル化します。トランスミッション システムは、トルク ドライブ、ドライブ シャフト、クラッチ、および出力シャフトに接続されている高速ギアと低速ギアで構成されます。

mdl = 'pdmTransmissionCasing';
open_system(mdl)

トランスミッション システムにはケーシングの振動を監視する振動センサーが含まれています。ケーシング モデルはシャフトの角変位をケーシングの線形変位に変換します。ケーシングは質量ばねダンパー システムとしてモデル化され、ケーシングからの振動 (ケーシング加速度) が測定されます。

open_system([mdl '/Casing'])

故障のモデル化

トランスミッション システムには、振動センサー ドリフト、ギア歯の故障、およびシャフト摩擦の各故障モデルが含まれています。センサー ドリフトの故障は、センサー モデルにオフセットを導入することで簡単にモデル化されます。オフセットはモデル変数 SDrift によって制御されます。SDrift=0 はセンサー故障がないことを意味する点に注意してください。

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

シャフト摩耗の故障はバリアント サブシステムによってモデル化されます。ここではサブシステム バリアントがシャフト減衰を変更しますが、バリアント サブシステムを使用してシャフト モデルの実装を完全に変えることもできます。選択したバリアントはモデル変数 ShaftWear によって制御されます。 ShaftWear=0 はシャフト故障がないことを意味する点に注意してください。

open_system([mdl '/Shaft'])

open_system([mdl,'/Shaft/Output Shaft'])

ギア歯の故障は、ドライブ シャフトの回転における固定位置に外乱トルクを挿入することでモデル化されます。シャフト位置はラジアン単位で測定され、シャフト位置が 0 近傍の小さいウィンドウ内にあるとシャフトに外乱力が適用されます。外乱の大きさはモデル変数 ToothFaultGain によって制御されます。ToothFaultGain=0 はギア歯の故障がないことを意味する点に注意してください。

故障データと健全データのシミュレーション

トランスミッション モデルは、センサー ドリフト、ギア歯、およびシャフト摩耗という 3 つの異なる故障タイプの存在と重大度を制御する変数を使って構成されます。モデル変数 SDriftToothFaultGain、および 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

同様に、各モデル変数の乱数値の組み合わせを作成します。3 つの故障タイプのうち一部のみを表す組み合わせができるように、必ず 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 は、ログ データをファイルに保存し、信号のログに timetable 形式を使用して、保存ファイルに Simulink.SimulationInput オブジェクトを格納するようモデルを構成します。関数 generateSimulationEnsemble は、シミュレーションが正常に完了したかどうかを示すステータス フラグを返します。

上記のコードは、グリッドの変数値から 110 のシミュレーション入力、ランダム変数の値から 98 のシミュレーション入力を作成し、合計 208 のシミュレーションが提供されています。これら 208 のシミュレーションを並列で実行するには標準的デスクトップで約 10 分を要し、約 10 GB のデータが生成されます。便宜上、最初の 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 信号には、ドライブおよび負荷シャフトの回転のパルスが含まれますが、解析 (具体的には、時間同期平均化) にはシャフト回転の時間が必要です。次のコードは、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 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"];

上記のコードはアンサンブルの 1 つのメンバーを処理するために使用しました。アンサンブルのすべてのメンバーを処理するには、上のコードを関数 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.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')

ピーク振幅に対応する周波数は、さまざまな故障タイプの分類に役立つ特徴となる可能性があります。下記のコードは、すべてのアンサンブル メンバーについて前述の特徴を計算します (この解析の実行には標準的デスクトップで約 30 分かかります。アンサンブルの 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)

すべてのアンサンブル メンバーの特徴データを 1 つの table に集めます。

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    
      ⋮

センサー ドリフト故障について考えます。上で計算したすべての特徴を予測子とし、センサー ドリフト故障のラベル (true または false の値) を応答として、fscnca コマンドを使用します。fscnca コマンドは各特徴の重みを返しますが、重みの大きい特徴は、応答を予測する上での重要度が高くなります。センサー ドリフト故障の場合、重みは、2 つの特徴が重要な予測子となり (信号の累積和範囲とスペクトル尖度のピーク周波数)、残りの特徴はほとんど影響しないことを示しています。

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 が得られます。特徴テーブルを、学習に使用されるメンバーとテストおよび検証に使用されるメンバーに分割します。学習用メンバーを使用して、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');   

混同行列は、分類器によって故障なしの条件はすべて正しく分類されるものの、故障と予期される条件の 1 つが故障なしとして誤分類されることを示しています。分類器で使用される特徴の数を増やすことにより、性能をさらに改善することができます。

まとめ

この例では、シミュレーション アンサンブルを使ってシミュレーション データのクリーン アップと特徴の抽出を行い、Simulink から故障データを生成するためのワークフローを説明しました。その後、抽出した特徴を使用して、異なった故障タイプ用の分類器を作成しました。

参考

関連するトピック