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 つの異なる故障タイプの存在と重大度を制御する変数を使って構成されます。モデル変数 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

同様に、各モデル変数の乱数値の組み合わせを作成します。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 のシミュレーションを並列実行すると、標準的なデスクトップでは 2 時間以上かかる場合があり、約 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

[21-Jul-2022 15:36:27] Checking for availability of parallel pool...
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 6).
[21-Jul-2022 15:37:34] Starting Simulink on parallel workers...
[21-Jul-2022 15:38:00] Configuring simulation cache folder on parallel workers...
[21-Jul-2022 15:38:01] Transferring base workspace variables used in the model to parallel workers...
[21-Jul-2022 15:38:01] Loading model on parallel workers...
[21-Jul-2022 15:38:20] Running simulations...
[21-Jul-2022 15:38:55] Completed 1 of 208 simulation runs
[21-Jul-2022 15:38:56] Completed 2 of 208 simulation runs
[21-Jul-2022 15:38:57] Completed 3 of 208 simulation runs
[21-Jul-2022 15:38:58] Completed 4 of 208 simulation runs
[21-Jul-2022 15:38:59] Completed 5 of 208 simulation runs
[21-Jul-2022 15:39:00] Completed 6 of 208 simulation runs
[21-Jul-2022 15:39:09] Completed 7 of 208 simulation runs
[21-Jul-2022 15:39:09] Completed 8 of 208 simulation runs
[21-Jul-2022 15:39:10] Completed 9 of 208 simulation runs
[21-Jul-2022 15:39:11] Completed 10 of 208 simulation runs
[21-Jul-2022 15:39:12] Completed 11 of 208 simulation runs
[21-Jul-2022 15:39:13] Completed 12 of 208 simulation runs
[21-Jul-2022 15:39:21] Completed 13 of 208 simulation runs
[21-Jul-2022 15:39:22] Completed 14 of 208 simulation runs
[21-Jul-2022 15:39:23] Completed 15 of 208 simulation runs
[21-Jul-2022 15:39:23] Completed 16 of 208 simulation runs
[21-Jul-2022 15:39:24] Completed 17 of 208 simulation runs
[21-Jul-2022 15:39:25] Completed 18 of 208 simulation runs
[21-Jul-2022 15:39:34] Completed 19 of 208 simulation runs
[21-Jul-2022 15:39:35] Completed 20 of 208 simulation runs
[21-Jul-2022 15:39:35] Completed 21 of 208 simulation runs
[21-Jul-2022 15:39:36] Completed 22 of 208 simulation runs
[21-Jul-2022 15:39:37] Completed 23 of 208 simulation runs
[21-Jul-2022 15:39:38] Completed 24 of 208 simulation runs
[21-Jul-2022 15:39:47] Completed 25 of 208 simulation runs
[21-Jul-2022 15:39:47] Completed 26 of 208 simulation runs
[21-Jul-2022 15:39:48] Completed 27 of 208 simulation runs
[21-Jul-2022 15:39:49] Completed 28 of 208 simulation runs
[21-Jul-2022 15:39:50] Completed 29 of 208 simulation runs
[21-Jul-2022 15:39:51] Completed 30 of 208 simulation runs
[21-Jul-2022 15:40:00] Completed 31 of 208 simulation runs
[21-Jul-2022 15:40:00] Completed 32 of 208 simulation runs
[21-Jul-2022 15:40:01] Completed 33 of 208 simulation runs
[21-Jul-2022 15:40:02] Completed 34 of 208 simulation runs
[21-Jul-2022 15:40:03] Completed 35 of 208 simulation runs
[21-Jul-2022 15:40:05] Completed 36 of 208 simulation runs
[21-Jul-2022 15:40:13] Completed 37 of 208 simulation runs
[21-Jul-2022 15:40:14] Completed 38 of 208 simulation runs
[21-Jul-2022 15:40:15] Completed 39 of 208 simulation runs
[21-Jul-2022 15:40:16] Completed 40 of 208 simulation runs
[21-Jul-2022 15:40:17] Completed 41 of 208 simulation runs
[21-Jul-2022 15:40:18] Completed 42 of 208 simulation runs
[21-Jul-2022 15:40:26] Completed 43 of 208 simulation runs
[21-Jul-2022 15:40:27] Completed 44 of 208 simulation runs
[21-Jul-2022 15:40:27] Completed 45 of 208 simulation runs
[21-Jul-2022 15:40:28] Completed 46 of 208 simulation runs
[21-Jul-2022 15:40:29] Completed 47 of 208 simulation runs
[21-Jul-2022 15:40:30] Completed 48 of 208 simulation runs
[21-Jul-2022 15:40:38] Completed 49 of 208 simulation runs
[21-Jul-2022 15:40:39] Completed 50 of 208 simulation runs
[21-Jul-2022 15:40:40] Completed 51 of 208 simulation runs
[21-Jul-2022 15:40:41] Completed 52 of 208 simulation runs
[21-Jul-2022 15:40:42] Completed 53 of 208 simulation runs
[21-Jul-2022 15:40:43] Completed 54 of 208 simulation runs
[21-Jul-2022 15:40:52] Completed 55 of 208 simulation runs
[21-Jul-2022 15:40:55] Completed 56 of 208 simulation runs
[21-Jul-2022 15:40:56] Completed 57 of 208 simulation runs
[21-Jul-2022 15:40:57] Completed 58 of 208 simulation runs
[21-Jul-2022 15:40:58] Completed 59 of 208 simulation runs
[21-Jul-2022 15:40:59] Completed 60 of 208 simulation runs
[21-Jul-2022 15:41:09] Completed 61 of 208 simulation runs
[21-Jul-2022 15:41:10] Completed 62 of 208 simulation runs
[21-Jul-2022 15:41:11] Completed 63 of 208 simulation runs
[21-Jul-2022 15:41:12] Completed 64 of 208 simulation runs
[21-Jul-2022 15:41:13] Completed 65 of 208 simulation runs
[21-Jul-2022 15:41:14] Completed 66 of 208 simulation runs
[21-Jul-2022 15:41:20] Completed 67 of 208 simulation runs
[21-Jul-2022 15:41:21] Completed 68 of 208 simulation runs
[21-Jul-2022 15:41:22] Completed 69 of 208 simulation runs
[21-Jul-2022 15:41:22] Completed 70 of 208 simulation runs
[21-Jul-2022 15:41:23] Completed 71 of 208 simulation runs
[21-Jul-2022 15:41:24] Completed 72 of 208 simulation runs
[21-Jul-2022 15:41:33] Completed 73 of 208 simulation runs
[21-Jul-2022 15:41:33] Completed 74 of 208 simulation runs
[21-Jul-2022 15:41:34] Completed 75 of 208 simulation runs
[21-Jul-2022 15:41:35] Completed 76 of 208 simulation runs
[21-Jul-2022 15:41:36] Completed 77 of 208 simulation runs
[21-Jul-2022 15:41:37] Completed 78 of 208 simulation runs
[21-Jul-2022 15:41:45] Completed 79 of 208 simulation runs
[21-Jul-2022 15:41:46] Completed 80 of 208 simulation runs
[21-Jul-2022 15:41:47] Completed 81 of 208 simulation runs
[21-Jul-2022 15:41:48] Completed 82 of 208 simulation runs
[21-Jul-2022 15:41:49] Completed 83 of 208 simulation runs
[21-Jul-2022 15:41:49] Completed 84 of 208 simulation runs
[21-Jul-2022 15:41:57] Completed 85 of 208 simulation runs
[21-Jul-2022 15:41:58] Completed 86 of 208 simulation runs
[21-Jul-2022 15:41:59] Completed 87 of 208 simulation runs
[21-Jul-2022 15:42:00] Completed 88 of 208 simulation runs
[21-Jul-2022 15:42:01] Completed 89 of 208 simulation runs
[21-Jul-2022 15:42:02] Completed 90 of 208 simulation runs
[21-Jul-2022 15:42:10] Completed 91 of 208 simulation runs
[21-Jul-2022 15:42:11] Completed 92 of 208 simulation runs
[21-Jul-2022 15:42:12] Completed 93 of 208 simulation runs
[21-Jul-2022 15:42:13] Completed 94 of 208 simulation runs
[21-Jul-2022 15:42:14] Completed 95 of 208 simulation runs
[21-Jul-2022 15:42:15] Completed 96 of 208 simulation runs
[21-Jul-2022 15:42:23] Completed 97 of 208 simulation runs
[21-Jul-2022 15:42:23] Completed 98 of 208 simulation runs
[21-Jul-2022 15:42:24] Completed 99 of 208 simulation runs
[21-Jul-2022 15:42:25] Completed 100 of 208 simulation runs
[21-Jul-2022 15:42:26] Completed 101 of 208 simulation runs
[21-Jul-2022 15:42:27] Completed 102 of 208 simulation runs
[21-Jul-2022 15:42:35] Completed 103 of 208 simulation runs
[21-Jul-2022 15:42:36] Completed 104 of 208 simulation runs
[21-Jul-2022 15:42:37] Completed 105 of 208 simulation runs
[21-Jul-2022 15:42:38] Completed 106 of 208 simulation runs
[21-Jul-2022 15:42:38] Completed 107 of 208 simulation runs
[21-Jul-2022 15:42:39] Completed 108 of 208 simulation runs
[21-Jul-2022 15:42:47] Completed 109 of 208 simulation runs
[21-Jul-2022 15:42:48] Completed 110 of 208 simulation runs
[21-Jul-2022 15:42:50] Completed 111 of 208 simulation runs
[21-Jul-2022 15:42:51] Completed 112 of 208 simulation runs
[21-Jul-2022 15:42:52] Completed 113 of 208 simulation runs
[21-Jul-2022 15:42:53] Completed 114 of 208 simulation runs
[21-Jul-2022 15:43:01] Completed 115 of 208 simulation runs
[21-Jul-2022 15:43:02] Completed 116 of 208 simulation runs
[21-Jul-2022 15:43:03] Completed 117 of 208 simulation runs
[21-Jul-2022 15:43:04] Completed 118 of 208 simulation runs
[21-Jul-2022 15:43:05] Completed 119 of 208 simulation runs
[21-Jul-2022 15:43:06] Completed 120 of 208 simulation runs
[21-Jul-2022 15:43:13] Completed 121 of 208 simulation runs
[21-Jul-2022 15:43:14] Completed 122 of 208 simulation runs
[21-Jul-2022 15:43:14] Completed 123 of 208 simulation runs
[21-Jul-2022 15:43:15] Completed 124 of 208 simulation runs
[21-Jul-2022 15:43:15] Completed 125 of 208 simulation runs
[21-Jul-2022 15:43:16] Completed 126 of 208 simulation runs
[21-Jul-2022 15:43:28] Completed 127 of 208 simulation runs
[21-Jul-2022 15:43:29] Completed 128 of 208 simulation runs
[21-Jul-2022 15:43:30] Completed 129 of 208 simulation runs
[21-Jul-2022 15:43:31] Completed 130 of 208 simulation runs
[21-Jul-2022 15:43:32] Completed 131 of 208 simulation runs
[21-Jul-2022 15:43:33] Completed 132 of 208 simulation runs
[21-Jul-2022 15:43:38] Completed 133 of 208 simulation runs
[21-Jul-2022 15:43:40] Completed 134 of 208 simulation runs
[21-Jul-2022 15:43:41] Completed 135 of 208 simulation runs
[21-Jul-2022 15:43:42] Completed 136 of 208 simulation runs
[21-Jul-2022 15:43:43] Completed 137 of 208 simulation runs
[21-Jul-2022 15:43:44] Completed 138 of 208 simulation runs
[21-Jul-2022 15:43:51] Completed 139 of 208 simulation runs
[21-Jul-2022 15:43:52] Completed 140 of 208 simulation runs
[21-Jul-2022 15:43:53] Completed 141 of 208 simulation runs
[21-Jul-2022 15:43:54] Completed 142 of 208 simulation runs
[21-Jul-2022 15:43:55] Completed 143 of 208 simulation runs
[21-Jul-2022 15:43:56] Completed 144 of 208 simulation runs
[21-Jul-2022 15:44:03] Completed 145 of 208 simulation runs
[21-Jul-2022 15:44:05] Completed 146 of 208 simulation runs
[21-Jul-2022 15:44:06] Completed 147 of 208 simulation runs
[21-Jul-2022 15:44:07] Completed 148 of 208 simulation runs
[21-Jul-2022 15:44:08] Completed 149 of 208 simulation runs
[21-Jul-2022 15:44:09] Completed 150 of 208 simulation runs
[21-Jul-2022 15:44:16] Completed 151 of 208 simulation runs
[21-Jul-2022 15:44:17] Completed 152 of 208 simulation runs
[21-Jul-2022 15:44:18] Completed 153 of 208 simulation runs
[21-Jul-2022 15:44:18] Completed 154 of 208 simulation runs
[21-Jul-2022 15:44:19] Completed 155 of 208 simulation runs
[21-Jul-2022 15:44:20] Completed 156 of 208 simulation runs
[21-Jul-2022 15:44:28] Completed 157 of 208 simulation runs
[21-Jul-2022 15:44:29] Completed 158 of 208 simulation runs
[21-Jul-2022 15:44:29] Completed 159 of 208 simulation runs
[21-Jul-2022 15:44:31] Completed 160 of 208 simulation runs
[21-Jul-2022 15:44:32] Completed 161 of 208 simulation runs
[21-Jul-2022 15:44:33] Completed 162 of 208 simulation runs
[21-Jul-2022 15:44:42] Completed 163 of 208 simulation runs
[21-Jul-2022 15:44:42] Completed 164 of 208 simulation runs
[21-Jul-2022 15:44:43] Completed 165 of 208 simulation runs
[21-Jul-2022 15:44:44] Completed 166 of 208 simulation runs
[21-Jul-2022 15:44:45] Completed 167 of 208 simulation runs
[21-Jul-2022 15:44:46] Completed 168 of 208 simulation runs
[21-Jul-2022 15:44:54] Completed 169 of 208 simulation runs
[21-Jul-2022 15:44:55] Completed 170 of 208 simulation runs
[21-Jul-2022 15:44:56] Completed 171 of 208 simulation runs
[21-Jul-2022 15:44:57] Completed 172 of 208 simulation runs
[21-Jul-2022 15:44:57] Completed 173 of 208 simulation runs
[21-Jul-2022 15:44:58] Completed 174 of 208 simulation runs
[21-Jul-2022 15:45:07] Completed 175 of 208 simulation runs
[21-Jul-2022 15:45:08] Completed 176 of 208 simulation runs
[21-Jul-2022 15:45:08] Completed 177 of 208 simulation runs
[21-Jul-2022 15:45:09] Completed 178 of 208 simulation runs
[21-Jul-2022 15:45:10] Completed 179 of 208 simulation runs
[21-Jul-2022 15:45:11] Completed 180 of 208 simulation runs
[21-Jul-2022 15:45:20] Completed 181 of 208 simulation runs
[21-Jul-2022 15:45:21] Completed 182 of 208 simulation runs
[21-Jul-2022 15:45:22] Completed 183 of 208 simulation runs
[21-Jul-2022 15:45:23] Completed 184 of 208 simulation runs
[21-Jul-2022 15:45:24] Completed 185 of 208 simulation runs
[21-Jul-2022 15:45:25] Completed 186 of 208 simulation runs
[21-Jul-2022 15:45:32] Completed 187 of 208 simulation runs
[21-Jul-2022 15:45:33] Completed 188 of 208 simulation runs
[21-Jul-2022 15:45:34] Completed 189 of 208 simulation runs
[21-Jul-2022 15:45:35] Completed 190 of 208 simulation runs
[21-Jul-2022 15:45:36] Completed 191 of 208 simulation runs
[21-Jul-2022 15:45:37] Completed 192 of 208 simulation runs
[21-Jul-2022 15:45:45] Completed 193 of 208 simulation runs
[21-Jul-2022 15:45:45] Completed 194 of 208 simulation runs
[21-Jul-2022 15:45:46] Completed 195 of 208 simulation runs
[21-Jul-2022 15:45:47] Completed 196 of 208 simulation runs
[21-Jul-2022 15:45:48] Completed 197 of 208 simulation runs
[21-Jul-2022 15:45:49] Completed 198 of 208 simulation runs
[21-Jul-2022 15:45:57] Completed 199 of 208 simulation runs
[21-Jul-2022 15:45:58] Completed 200 of 208 simulation runs
[21-Jul-2022 15:45:59] Completed 201 of 208 simulation runs
[21-Jul-2022 15:46:00] Completed 202 of 208 simulation runs
[21-Jul-2022 15:46:01] Completed 203 of 208 simulation runs
[21-Jul-2022 15:46:02] Completed 204 of 208 simulation runs
[21-Jul-2022 15:46:09] Completed 205 of 208 simulation runs
[21-Jul-2022 15:46:10] Completed 206 of 208 simulation runs
[21-Jul-2022 15:46:11] Completed 207 of 208 simulation runs
[21-Jul-2022 15:46:12] Completed 208 of 208 simulation runs
[21-Jul-2022 15:46:12] 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]
                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"];

上記のコードはアンサンブルの 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.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)

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

featureData = gather(tall(ens))

センサードリフト故障について考えます。上で計算したすべての特徴を予測子とし、センサードリフト故障のラベル (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
idxSelectedFeature = featureAnalysis.FeatureWeights > 0.1;
classifySD = [featureData(:,idxSelectedFeature), featureData(:,idxResponse)]

累積和範囲のグループ化されたヒストグラムから、この特徴がなぜセンサードリフト故障の重要な予測子であるかを把握することができます。

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)]

信号のリアプノフ指数のグループ化されたヒストグラムには、この特徴だけでは適切な予測子とならない理由が示されます。

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