このページの翻訳は最新ではありません。ここをクリックして、英語の最新版を参照してください。
シミュレーション データを使用した複数クラス故障検出
この例では、Simulink モデルを使用して健全な状態のデータと故障状態のデータを生成する方法を説明します。このデータを使用して、さまざまな故障の組み合わせを検出する複数クラス分類器を開発します。この例は 3 重往復ポンプ モデルを使用し、漏れ、閉塞、およびベアリングの故障を含みます。
モデルの設定
この例では、zip ファイルに保存されている多くのサポート ファイルを使用します。ファイルを解凍してサポート ファイルにアクセスし、モデル パラメーターを読み込み、往復ポンプ ライブラリを作成します。
if ~exist('+mech_hydro_forcesPS','dir') unzip('pdmRecipPump_supportingfiles.zip') end % Load Parameters pdmRecipPump_Parameters %Pump CAT_Pump_1051_DataFile_imported %CAD % Create Simscape library if needed if exist('mech_hydro_forcesPS_Lib','file')~=4 ssc_build mech_hydro_forcesPS end
往復ポンプ モデル
往復ポンプは、電気モーター、ポンプ ハウジング、ポンプ クランク、およびポンプ プランジャーで構成されています。
mdl = 'pdmRecipPump';
open_system(mdl)
open_system([mdl,'/Pump'])
ポンプ モデルは、シリンダーの漏れ、吸込口の閉塞、およびベアリング摩擦の増加という 3 タイプの故障をモデル化するように設定されます。これらの故障はワークスペース変数としてパラメーター化され、ポンプ ブロック ダイアログによって設定されます。
故障データと健全データのシミュレーション
3 つの故障タイプそれぞれについて、故障なしから重大な故障までの故障重大度を表す値の配列を作成します。
% Define fault parameter variations numParValues = 10; leak_area_set_factor = linspace(0.00,0.036,numParValues); leak_area_set = leak_area_set_factor*TRP_Par.Check_Valve.In.Max_Area; leak_area_set = max(leak_area_set,1e-9); % Leakage area cannot be 0 blockinfactor_set = linspace(0.8,0.53,numParValues); bearingfactor_set = linspace(0,6e-4,numParValues);
ポンプ モデルはノイズを含めるように設定されているため、同じ故障パラメーター値でモデルを実行すると、結果としてさまざまなシミュレーション出力が得られます。これは、同じ故障状態と重大度に対して複数のシミュレーション結果が可能であることを意味するので、分類器の開発に役立ちます。そのような結果が得られるシミュレーションを設定するには、故障パラメーター値のベクトルを作成します。ここで値は、故障なし、1 つの故障、2 つの故障の組み合わせ、および 3 つの故障の組み合わせを表します。故障なし、1 つの故障、などの各グループについて、上記で定義されている故障パラメーター値から、故障値の 125 の組み合わせを作成します。これにより故障パラメーター値の合計 1000 の組み合わせが得られます。これら 1000 のシミュレーションを並列で実行する場合、標準のデスクトップでは約 1 時間かかり、約 620 MB のデータが生成されることに注意してください。シミュレーション時間を短縮させるには、runAll = true を runAll = false に変更して、故障の組み合わせの数を 20 に減らします。データセットが大きいほど、よりロバストな分類器となることに注意してください。
% Set number of elements in each fault group runAll = false; if runAll % Create a large dataset to build a robust classifier nPerGroup = 125; else % Create a smaller dataset to reduce simulation time nPerGroup = 20; %#ok<UNRCH> end % No fault simulations leakArea = repmat(leak_area_set(1),nPerGroup,1); blockingFactor = repmat(blockinfactor_set(1),nPerGroup,1); bearingFactor = repmat(bearingfactor_set(1),nPerGroup,1); % Single fault simulations idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idx)']; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;blockinfactor_set(idx)']; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idx = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idx)']; % Double fault simulations idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;blockinfactor_set(idxB)']; bearingFactor = [bearingFactor;repmat(bearingfactor_set(1),nPerGroup,1)]; idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;repmat(blockinfactor_set(1),nPerGroup,1)]; bearingFactor = [bearingFactor;bearingfactor_set(idxB)']; idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; repmat(leak_area_set(1),nPerGroup,1)]; blockingFactor = [blockingFactor;blockinfactor_set(idxA)']; bearingFactor = [bearingFactor;bearingfactor_set(idxB)']; % Triple fault simulations idxA = ceil(10*rand(nPerGroup,1)); idxB = ceil(10*rand(nPerGroup,1)); idxC = ceil(10*rand(nPerGroup,1)); leakArea = [leakArea; leak_area_set(idxA)']; blockingFactor = [blockingFactor;blockinfactor_set(idxB)']; bearingFactor = [bearingFactor;bearingfactor_set(idxC)'];
故障パラメーターの組み合わせを使用して Simulink.SimulationInput オブジェクトを作成します。各シミュレーション入力につき、異なる結果が生成されるようにランダム シードを必ず異なる値に設定します。
for ct = numel(leakArea):-1:1 simInput(ct) = Simulink.SimulationInput(mdl); simInput(ct) = setVariable(simInput(ct),'leak_cyl_area_WKSP',leakArea(ct)); simInput(ct) = setVariable(simInput(ct),'block_in_factor_WKSP',blockingFactor(ct)); simInput(ct) = setVariable(simInput(ct),'bearing_fault_frict_WKSP',bearingFactor(ct)); simInput(ct) = setVariable(simInput(ct),'noise_seed_offset_WKSP',ct-1); end
関数 generateSimulationEnsemble を使用して、上で定義された Simulink.SimulationInput オブジェクトによって定義されるシミュレーションを実行し、ローカル サブフォルダーに結果を格納します。次に、格納された結果から simulationEnsembleDatastore を作成します。
% Run the simulation and create an ensemble to manage the simulation % results if isfolder('./Data') % Delete existing mat files delete('./Data/*.mat') end [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);
[08-Sep-2020 21:01:46] Checking for availability of parallel pool... Starting parallel pool (parpool) using the 'local' profile ... Connected to the parallel pool (number of workers: 6). [08-Sep-2020 21:02:27] Starting Simulink on parallel workers... [08-Sep-2020 21:02:55] Configuring simulation cache folder on parallel workers... [08-Sep-2020 21:02:55] Loading model on parallel workers... [08-Sep-2020 21:03:02] Transferring base workspace variables used in the model to parallel workers... [08-Sep-2020 21:03:08] Running simulations... [08-Sep-2020 21:04:21] Completed 1 of 160 simulation runs [08-Sep-2020 21:04:21] Completed 2 of 160 simulation runs [08-Sep-2020 21:04:22] Completed 3 of 160 simulation runs [08-Sep-2020 21:04:23] Completed 4 of 160 simulation runs [08-Sep-2020 21:04:24] Completed 5 of 160 simulation runs [08-Sep-2020 21:04:24] Completed 6 of 160 simulation runs [08-Sep-2020 21:04:47] Completed 7 of 160 simulation runs [08-Sep-2020 21:04:48] Completed 8 of 160 simulation runs [08-Sep-2020 21:04:49] Completed 9 of 160 simulation runs [08-Sep-2020 21:04:50] Completed 10 of 160 simulation runs [08-Sep-2020 21:04:52] Completed 11 of 160 simulation runs [08-Sep-2020 21:04:53] Completed 12 of 160 simulation runs [08-Sep-2020 21:05:14] Completed 13 of 160 simulation runs [08-Sep-2020 21:05:15] Completed 14 of 160 simulation runs [08-Sep-2020 21:05:17] Completed 15 of 160 simulation runs [08-Sep-2020 21:05:18] Completed 16 of 160 simulation runs [08-Sep-2020 21:05:20] Completed 17 of 160 simulation runs [08-Sep-2020 21:05:21] Completed 18 of 160 simulation runs [08-Sep-2020 21:05:41] Completed 19 of 160 simulation runs [08-Sep-2020 21:05:43] Completed 20 of 160 simulation runs [08-Sep-2020 21:05:45] Completed 21 of 160 simulation runs [08-Sep-2020 21:05:47] Completed 22 of 160 simulation runs [08-Sep-2020 21:05:49] Completed 23 of 160 simulation runs [08-Sep-2020 21:05:51] Completed 24 of 160 simulation runs [08-Sep-2020 21:06:08] Completed 25 of 160 simulation runs [08-Sep-2020 21:06:10] Completed 26 of 160 simulation runs [08-Sep-2020 21:06:13] Completed 27 of 160 simulation runs [08-Sep-2020 21:06:16] Completed 28 of 160 simulation runs [08-Sep-2020 21:06:18] Completed 29 of 160 simulation runs [08-Sep-2020 21:06:21] Completed 30 of 160 simulation runs [08-Sep-2020 21:06:34] Completed 31 of 160 simulation runs [08-Sep-2020 21:06:38] Completed 32 of 160 simulation runs [08-Sep-2020 21:06:41] Completed 33 of 160 simulation runs [08-Sep-2020 21:06:44] Completed 34 of 160 simulation runs [08-Sep-2020 21:06:48] Completed 35 of 160 simulation runs [08-Sep-2020 21:06:51] Completed 36 of 160 simulation runs [08-Sep-2020 21:07:01] Completed 37 of 160 simulation runs [08-Sep-2020 21:07:05] Completed 38 of 160 simulation runs [08-Sep-2020 21:07:09] Completed 39 of 160 simulation runs [08-Sep-2020 21:07:13] Completed 40 of 160 simulation runs [08-Sep-2020 21:07:17] Completed 41 of 160 simulation runs [08-Sep-2020 21:07:21] Completed 42 of 160 simulation runs [08-Sep-2020 21:07:32] Completed 43 of 160 simulation runs [08-Sep-2020 21:07:37] Completed 44 of 160 simulation runs [08-Sep-2020 21:07:41] Completed 45 of 160 simulation runs [08-Sep-2020 21:07:46] Completed 46 of 160 simulation runs [08-Sep-2020 21:07:50] Completed 47 of 160 simulation runs [08-Sep-2020 21:07:58] Completed 48 of 160 simulation runs [08-Sep-2020 21:08:04] Completed 49 of 160 simulation runs [08-Sep-2020 21:08:10] Completed 50 of 160 simulation runs [08-Sep-2020 21:08:16] Completed 51 of 160 simulation runs [08-Sep-2020 21:08:23] Completed 52 of 160 simulation runs [08-Sep-2020 21:08:33] Completed 53 of 160 simulation runs [08-Sep-2020 21:08:41] Completed 54 of 160 simulation runs [08-Sep-2020 21:08:51] Completed 55 of 160 simulation runs [08-Sep-2020 21:09:01] Completed 56 of 160 simulation runs [08-Sep-2020 21:09:11] Completed 57 of 160 simulation runs [08-Sep-2020 21:09:22] Completed 58 of 160 simulation runs [08-Sep-2020 21:09:34] Completed 59 of 160 simulation runs [08-Sep-2020 21:09:45] Completed 60 of 160 simulation runs [08-Sep-2020 21:09:57] Completed 61 of 160 simulation runs [08-Sep-2020 21:10:08] Completed 62 of 160 simulation runs [08-Sep-2020 21:10:16] Completed 63 of 160 simulation runs [08-Sep-2020 21:10:25] Completed 64 of 160 simulation runs [08-Sep-2020 21:10:36] Completed 65 of 160 simulation runs [08-Sep-2020 21:10:44] Completed 66 of 160 simulation runs [08-Sep-2020 21:10:57] Completed 67 of 160 simulation runs [08-Sep-2020 21:11:10] Completed 68 of 160 simulation runs [08-Sep-2020 21:11:21] Completed 69 of 160 simulation runs [08-Sep-2020 21:11:31] Completed 70 of 160 simulation runs [08-Sep-2020 21:11:44] Completed 71 of 160 simulation runs [08-Sep-2020 21:11:55] Completed 72 of 160 simulation runs [08-Sep-2020 21:12:11] Completed 73 of 160 simulation runs [08-Sep-2020 21:12:28] Completed 74 of 160 simulation runs [08-Sep-2020 21:12:43] Completed 75 of 160 simulation runs [08-Sep-2020 21:12:59] Completed 76 of 160 simulation runs [08-Sep-2020 21:13:12] Completed 77 of 160 simulation runs [08-Sep-2020 21:13:27] Completed 78 of 160 simulation runs [08-Sep-2020 21:13:43] Completed 79 of 160 simulation runs [08-Sep-2020 21:13:58] Completed 80 of 160 simulation runs [08-Sep-2020 21:14:18] Completed 81 of 160 simulation runs [08-Sep-2020 21:14:30] Completed 82 of 160 simulation runs [08-Sep-2020 21:14:36] Completed 83 of 160 simulation runs [08-Sep-2020 21:14:42] Completed 84 of 160 simulation runs [08-Sep-2020 21:14:47] Completed 85 of 160 simulation runs [08-Sep-2020 21:14:53] Completed 86 of 160 simulation runs [08-Sep-2020 21:14:59] Completed 87 of 160 simulation runs [08-Sep-2020 21:15:05] Completed 88 of 160 simulation runs [08-Sep-2020 21:15:12] Completed 89 of 160 simulation runs [08-Sep-2020 21:15:18] Completed 90 of 160 simulation runs [08-Sep-2020 21:15:26] Completed 91 of 160 simulation runs [08-Sep-2020 21:15:33] Completed 92 of 160 simulation runs [08-Sep-2020 21:15:40] Completed 93 of 160 simulation runs [08-Sep-2020 21:15:46] Completed 94 of 160 simulation runs [08-Sep-2020 21:15:53] Completed 95 of 160 simulation runs [08-Sep-2020 21:16:00] Completed 96 of 160 simulation runs [08-Sep-2020 21:16:07] Completed 97 of 160 simulation runs [08-Sep-2020 21:16:14] Completed 98 of 160 simulation runs [08-Sep-2020 21:16:21] Completed 99 of 160 simulation runs [08-Sep-2020 21:16:29] Completed 100 of 160 simulation runs [08-Sep-2020 21:16:36] Completed 101 of 160 simulation runs [08-Sep-2020 21:16:41] Completed 102 of 160 simulation runs [08-Sep-2020 21:16:48] Completed 103 of 160 simulation runs [08-Sep-2020 21:16:53] Completed 104 of 160 simulation runs [08-Sep-2020 21:16:58] Completed 105 of 160 simulation runs [08-Sep-2020 21:17:02] Completed 106 of 160 simulation runs [08-Sep-2020 21:17:10] Completed 107 of 160 simulation runs [08-Sep-2020 21:17:15] Completed 108 of 160 simulation runs [08-Sep-2020 21:17:21] Completed 109 of 160 simulation runs [08-Sep-2020 21:17:26] Completed 110 of 160 simulation runs [08-Sep-2020 21:17:32] Completed 111 of 160 simulation runs [08-Sep-2020 21:17:37] Completed 112 of 160 simulation runs [08-Sep-2020 21:17:43] Completed 113 of 160 simulation runs [08-Sep-2020 21:17:48] Completed 114 of 160 simulation runs [08-Sep-2020 21:17:54] Completed 115 of 160 simulation runs [08-Sep-2020 21:18:00] Completed 116 of 160 simulation runs [08-Sep-2020 21:18:05] Completed 117 of 160 simulation runs [08-Sep-2020 21:18:12] Completed 118 of 160 simulation runs [08-Sep-2020 21:18:18] Completed 119 of 160 simulation runs [08-Sep-2020 21:18:24] Completed 120 of 160 simulation runs [08-Sep-2020 21:18:30] Completed 121 of 160 simulation runs [08-Sep-2020 21:18:36] Completed 122 of 160 simulation runs [08-Sep-2020 21:18:42] Completed 123 of 160 simulation runs [08-Sep-2020 21:18:49] Completed 124 of 160 simulation runs [08-Sep-2020 21:18:57] Completed 125 of 160 simulation runs [08-Sep-2020 21:19:03] Completed 126 of 160 simulation runs [08-Sep-2020 21:19:09] Completed 127 of 160 simulation runs [08-Sep-2020 21:19:14] Completed 128 of 160 simulation runs [08-Sep-2020 21:19:20] Completed 129 of 160 simulation runs [08-Sep-2020 21:19:28] Completed 130 of 160 simulation runs [08-Sep-2020 21:19:35] Completed 131 of 160 simulation runs [08-Sep-2020 21:19:41] Completed 132 of 160 simulation runs [08-Sep-2020 21:19:47] Completed 133 of 160 simulation runs [08-Sep-2020 21:19:54] Completed 134 of 160 simulation runs [08-Sep-2020 21:20:03] Completed 135 of 160 simulation runs [08-Sep-2020 21:20:12] Completed 136 of 160 simulation runs [08-Sep-2020 21:20:19] Completed 137 of 160 simulation runs [08-Sep-2020 21:20:27] Completed 138 of 160 simulation runs [08-Sep-2020 21:20:35] Completed 139 of 160 simulation runs [08-Sep-2020 21:20:43] Completed 140 of 160 simulation runs [08-Sep-2020 21:20:50] Completed 141 of 160 simulation runs [08-Sep-2020 21:20:56] Completed 142 of 160 simulation runs [08-Sep-2020 21:21:03] Completed 143 of 160 simulation runs [08-Sep-2020 21:21:13] Completed 144 of 160 simulation runs [08-Sep-2020 21:21:19] Completed 145 of 160 simulation runs [08-Sep-2020 21:21:25] Completed 146 of 160 simulation runs [08-Sep-2020 21:21:32] Completed 147 of 160 simulation runs [08-Sep-2020 21:21:42] Completed 148 of 160 simulation runs [08-Sep-2020 21:21:49] Completed 149 of 160 simulation runs [08-Sep-2020 21:21:57] Completed 150 of 160 simulation runs [08-Sep-2020 21:22:05] Completed 151 of 160 simulation runs [08-Sep-2020 21:22:14] Completed 152 of 160 simulation runs [08-Sep-2020 21:22:22] Completed 153 of 160 simulation runs [08-Sep-2020 21:22:31] Completed 154 of 160 simulation runs [08-Sep-2020 21:22:39] Completed 155 of 160 simulation runs [08-Sep-2020 21:22:47] Completed 156 of 160 simulation runs [08-Sep-2020 21:22:55] Completed 157 of 160 simulation runs [08-Sep-2020 21:23:03] Completed 158 of 160 simulation runs [08-Sep-2020 21:23:12] Completed 159 of 160 simulation runs [08-Sep-2020 21:23:21] Completed 160 of 160 simulation runs [08-Sep-2020 21:23:23] Cleaning up parallel workers...
ens = simulationEnsembleDatastore(fullfile('.','Data'));
シミュレーション結果からの特徴の処理と抽出
このモデルはポンプの吐出圧力、出流量、モーター速度、およびモーター電流をログに記録するよう設定されます。
ens.DataVariables
ans = 8×1 string
"SimulationInput"
"SimulationMetadata"
"iMotor_meas"
"pIn_meas"
"pOut_meas"
"qIn_meas"
"qOut_meas"
"wMotor_meas"
アンサンブルのメンバーごとに、ポンプ出流量を前処理し、ポンプ出流量に基づいて状態インジケーターを計算します。状態インジケーターは後で故障分類のために使用されます。前処理のため、出流量の最初の 0.8 秒を削除します。これにはシミュレーションおよびポンプ始動からの過渡状態が含まれるためです。前処理の一部として、出流量のパワー スペクトルを計算し、SimulationInput を使って故障変数の値を返します。
対象とする変数だけが読み取りで返されるようにアンサンブルを設定し、この例の最後で定義されている関数 preprocess を呼び出します。
ens.SelectedVariables = ["qOut_meas", "SimulationInput"]; reset(ens) data = read(ens)
data=1×2 table
qOut_meas SimulationInput
__________________ ______________________________
{2001×1 timetable} {1×1 Simulink.SimulationInput}
[flow,flowP,flowF,faultValues] = preprocess(data);
流量と流量スペクトルをプロットします。プロットされるデータは故障なしの状態のものです。
% Figure with nominal subplot(211); plot(flow.Time,flow.Data); subplot(212); semilogx(flowF,pow2db(flowP)); xlabel('Hz')
流量スペクトルは、期待された周波数での共振ピークを示しています。具体的には、ポンプのモーター速度は 950 rpm、つまり 15.833 Hz であり、ポンプはシリンダーを 3 つもつため、流量は 3*15.833 Hz つまり 47.5 Hz で基本波であり、47.5 Hz の倍数において高調波であると予想されます。流量スペクトルは、期待される共振ピークを明確に示しています。ポンプの 1 つのシリンダーに故障が生じると、ポンプのモーター速度 15.833 Hz およびその高調波で共振が発生します。
流量スペクトルおよび低速信号は、可能な状態インジケーターをある程度示唆しています。具体的には、平均や分散などの一般的な信号統計量、およびスペクトル特性です。ピーク振幅の周波数、15.833 Hz 周辺のエネルギー、47.5 Hz 周辺のエネルギー、および 100 Hz を超えるエネルギーなどの、期待される高調波に関連するスペクトル状態インジケーターが計算されます。スペクトル尖度ピークの周波数も計算されます。
状態インジケーターのデータ変数および故障変数値の状態変数を使用して、アンサンブルを設定します。その後、関数 extractCI を呼び出して特徴を計算し、writeToLastMemberRead コマンドを使用して特徴および故障変数の値をアンサンブルに追加します。関数 extractCI はこの例の最後で定義されています。
ens.DataVariables = [ens.DataVariables; ... "fPeak"; "pLow"; "pMid"; "pHigh"; "pKurtosis"; ... "qMean"; "qVar"; "qSkewness"; "qKurtosis"; ... "qPeak2Peak"; "qCrest"; "qRMS"; "qMAD"; "qCSRange"]; ens.ConditionVariables = ["LeakFault","BlockingFault","BearingFault"]; feat = extractCI(flow,flowP,flowF); dataToWrite = [faultValues, feat]; writeToLastMemberRead(ens,dataToWrite{:})
上記のコードは、シミュレーション アンサンブルの最初のメンバーの状態インジケーターを前処理して計算します。アンサンブルの hasdata コマンドを使用して、アンサンブルのすべてのメンバーについてこれを繰り返します。さまざまな故障状態下におけるシミュレーション結果を把握するには、アンサンブルを 100 要素ごとにプロットします。
%Figure with nominal and faults figure, subplot(211); lFlow = plot(flow.Time,flow.Data,'LineWidth',2); subplot(212); lFlowP = semilogx(flowF,pow2db(flowP),'LineWidth',2); xlabel('Hz') ct = 1; lColors = get(lFlow.Parent,'ColorOrder'); lIdx = 2; % Loop over all members in the ensemble, preprocess % and compute the features for each member while hasdata(ens) % Read member data data = read(ens); % Preprocess and extract features from the member data [flow,flowP,flowF,faultValues] = preprocess(data); feat = extractCI(flow,flowP,flowF); % Add the extracted feature values to the member data dataToWrite = [faultValues, feat]; writeToLastMemberRead(ens,dataToWrite{:}) % Plot member signal and spectrum for every 100th member if mod(ct,100) == 0 line('Parent',lFlow.Parent,'XData',flow.Time,'YData',flow.Data, ... 'Color', lColors(lIdx,:)); line('Parent',lFlowP.Parent,'XData',flowF,'YData',pow2db(flowP), ... 'Color', lColors(lIdx,:)); if lIdx == size(lColors,1) lIdx = 1; else lIdx = lIdx+1; end end ct = ct + 1; end
さまざまな故障状態および重大度の下において、期待される周波数でスペクトルに高調波が含まれることに注意してください。
ポンプの故障の検出と分類
前の節では、シミュレーション アンサンブルの全メンバーの流量信号から状態インジケーターを前処理して計算しました。これはさまざまな故障の組み合わせと重大度によるシミュレーション結果に対応します。状態インジケーターを使用して、ポンプの流量信号からポンプの故障を検出し分類することができます。
状態インジケーターを読み取るようにシミュレーション アンサンブルを設定し、tall コマンドと gather コマンドを使用してすべての状態インジケーターと故障変数値をメモリに読み込むことができます｡
% Get data to design a classifier. reset(ens) ens.SelectedVariables = [... "fPeak","pLow","pMid","pHigh","pKurtosis",... "qMean","qVar","qSkewness","qKurtosis",... "qPeak2Peak","qCrest","qRMS","qMAD","qCSRange",... "LeakFault","BlockingFault","BearingFault"]; idxLastFeature = 14; % Load the condition indicator data into memory data = gather(tall(ens));
Evaluating tall expression using the Parallel Pool 'local': - Pass 1 of 1: 0% complete Evaluation 0% complete
- Pass 1 of 1: Completed in 11 sec Evaluation completed in 11 sec
data(1:10,:)
ans=10×17 table
fPeak pLow pMid pHigh pKurtosis qMean qVar qSkewness qKurtosis qPeak2Peak qCrest qRMS qMAD qCSRange LeakFault BlockingFault BearingFault
______ _______ ______ ______ _________ ______ ______ _________ _________ __________ ______ ______ ______ ________ _________ _____________ ____________
43.909 0.84932 117.83 18.867 276.49 35.573 7.5235 -0.73064 2.778 13.86 1.1491 35.679 2.2319 42692 1e-09 0.8 0
43.909 0.46288 126.04 18.969 12.417 35.577 7.8766 -0.70939 2.63 13.336 1.1451 35.688 2.3235 42699 1e-09 0.8 0
44.03 66.853 151.07 23.624 230.13 34.357 9.4974 -0.54522 2.8046 15.243 1.1665 34.495 2.4703 41231 1.6e-06 0.71 0
40.699 6.0325 89.768 21.02 252.48 32.617 6.7424 -0.70295 2.7538 12.264 1.1457 32.721 2.1208 39139 4e-07 0.8 6.6667e-05
8.6012 14.366 29.295 9.8881 485.93 18.623 9.325 -0.02553 2.0546 14.315 1.3422 18.871 2.6236 22347 2.8e-06 0.8 0.0006
43.969 33.122 141.94 24.284 154.8 34.682 8.6405 -0.53746 2.7521 15.692 1.1955 34.807 2.3614 41620 1.2e-06 0.8 0
9.5096 25.219 31.388 13.27 65.397 21.352 6.7393 -0.25342 2.4885 13.704 1.2738 21.509 2.1357 25619 2e-06 0.8 0.00046667
26.83 1.114 31.794 15.963 434.6 21.432 3.6724 -0.44013 2.6938 10.304 1.2108 21.517 1.5589 25716 4e-07 0.8 0.00053333
9.0252 16.051 23.779 12.955 302.15 20.006 7.5907 -0.13753 2.2392 14.006 1.2837 20.195 2.3023 24010 2.4e-06 0.8 0.00053333
43.969 12.291 124.8 22.447 478.48 34.986 8.0954 -0.71245 2.9268 15.015 1.1585 35.101 2.2984 41987 8e-07 0.8 0
アンサンブル メンバーごとの故障変数値 (データ テーブルの 1 行) を故障フラグに変換し、これらの故障フラグを、各メンバーのさまざまな故障ステータスを捕捉する 1 つのフラグに組み合わせることが可能です。
% Convert the fault variable values to flags
data.LeakFlag = data.LeakFault > 1e-6;
data.BlockingFlag = data.BlockingFault < 0.8;
data.BearingFlag = data.BearingFault > 0;
data.CombinedFlag = data.LeakFlag+2*data.BlockingFlag+4*data.BearingFlag;状態インジケーターを入力として受け取り、組み合わせた故障フラグを返す分類器を作成します。2 次多項式カーネルを使用するサポート ベクター マシンを学習させます。cvpartition コマンドを使用して、アンサンブルのメンバーを学習用のセットと検証用のセットに分割します。
rng('default') % for reproducibility predictors = data(:,1:idxLastFeature); response = data.CombinedFlag; cvp = cvpartition(size(predictors,1),'KFold',5); % Create and train the classifier template = templateSVM(... 'KernelFunction', 'polynomial', ... 'PolynomialOrder', 2, ... 'KernelScale', 'auto', ... 'BoxConstraint', 1, ... 'Standardize', true); combinedClassifier = fitcecoc(... predictors(cvp.training(1),:), ... response(cvp.training(1),:), ... 'Learners', template, ... 'Coding', 'onevsone', ... 'ClassNames', [0; 1; 2; 3; 4; 5; 6; 7]);
検証データを使用して学習後の分類器の性能をチェックし、結果を混同プロットにプロットします。
% Check performance by computing and plotting the confusion matrix actualValue = response(cvp.test(1),:); predictedValue = predict(combinedClassifier, predictors(cvp.test(1),:)); confdata = confusionmat(actualValue,predictedValue); figure, labels = {'None', 'Leak','Blocking', 'Leak & Blocking', 'Bearing', ... 'Bearing & Leak', 'Bearing & Blocking', 'All'}; h = heatmap(confdata, ... 'YLabel', 'Actual leak fault', ... 'YDisplayLabels', labels, ... 'XLabel', 'Predicted fault', ... 'XDisplayLabels', labels, ... 'ColorbarVisible','off');
混同プロットは、故障の各組み合わせにつき、故障の組み合わせが正しく予測された回数 (プロットの対角エントリ) と、故障の組み合わせが誤って予測された回数 (非対角エントリ) を示しています。
混同プロットから、分類器は一部の故障状態を正しく分類しなかった (非対角項) ことがわかります。しかし、故障なしの状態は正しく予測されています。故障があっても故障なしの状態が予測された箇所 (1 列目) がいくつかありますが、それ以外には、厳密には正しい故障状態ではないにしても、故障が予測されています。全体的な検証の正確性は 84% で、故障があることを予測する正確性は 98% です。
% Compute the overall accuracy of the classifier
sum(diag(confdata))/sum(confdata(:))ans = 0.5000
% Compute the accuracy of the classifier at predicting % that there is a fault 1-sum(confdata(2:end,1))/sum(confdata(:))
ans = 0.9375
故障なしと予測され、故障が存在したケースについて調べます。まず、検証データで、実際の故障が閉塞故障でありながら、故障なしと予測されたケースを見つけます。
vData = data(cvp.test(1),:); b1 = (actualValue==2) & (predictedValue==0); fData = vData(b1,15:17)
fData=2×3 table
LeakFault BlockingFault BearingFault
_________ _____________ ____________
1e-09 0.74 0
8e-07 0.59 0
検証データで、実際の故障がベアリング故障でありながら、故障なしとして予測されたケースを見つけます。
b2 = (actualValue==4) & (predictedValue==0); vData(b2,15:17)
ans = 0×3 empty table
故障なしと予測されたにもかかわらず故障が存在したケースを調べてみると、これは閉塞故障値の 0.77 がそのノミナル値 0.8 に近い場合、またはベアリング故障値の 6.6e-5 がそのノミナル値 0 に近い場合に生じることが明らかになります。閉塞故障値が小さいケースのスペクトルをプロットし、故障なしの状態と比較すると、スペクトルがよく似ているために検出が困難になることがわかります。分類器を再学習させて、閉塞値の 0.77 を故障なしの状態として含めると、故障検出器の性能は大幅に改善されます。あるいは、追加のポンプ測定値を使用してさらに多くの情報を提供することにより、小さな閉塞故障を検出できるように機能を改善することができます。
% Configure the ensemble to only read the flow and fault variable values ens.SelectedVariables = ["qOut_meas","LeakFault","BlockingFault","BearingFault"]; reset(ens) % Load the ensemble member data into memory data = gather(tall(ens));
Evaluating tall expression using the Parallel Pool 'local': - Pass 1 of 1: 0% complete Evaluation 0% complete
- Pass 1 of 1: Completed in 9.7 sec Evaluation completed in 9.8 sec
% Look for the member that was incorrectly predicted, and % compute its power spectrum idx = ... data.BlockingFault == fData.BlockingFault(1) & ... data.LeakFault == fData.LeakFault(1) & ... data.BearingFault == fData.BearingFault(1); flow1 = data(idx,1); flow1 = flow1.qOut_meas{1}; [flow1P,flow1F] = pspectrum(flow1); % Look for a member that does not have any faults idx = ... data.BlockingFault == 0.8 & ... data.LeakFault == 1e-9 & ... data.BearingFault == 0; flow2 = data(idx,1); flow2 = flow2.qOut_meas{1}; [flow2P,flow2F] = pspectrum(flow2); % Plot the power spectra semilogx(... flow1F,pow2db(flow1P),... flow2F,pow2db(flow2P)); xlabel('Hz') legend('Small blocking fault','No fault')
まとめ
この例では Simulink モデルを使用して、往復ポンプの故障をモデル化し、さまざまな故障の組み合わせと重大度の下でモデルをシミュレートして、ポンプの出流量から状態インジケーターを抽出し、その状態インジケーターを使用してポンプの故障を検出するための分類器を学習させる方法を示しました。ここでは分類器を使用した故障検出の性能を調べ、小さな閉塞故障は故障なしの状態によく似ているため正確に検出できないことを確認しました。
サポート関数
function [flow,flowSpectrum,flowFrequencies,faultValues] = preprocess(data) % Helper function to preprocess the logged reciprocating pump data. % Remove the 1st 0.8 seconds of the flow signal tMin = seconds(0.8); flow = data.qOut_meas{1}; flow = flow(flow.Time >= tMin,:); flow.Time = flow.Time - flow.Time(1); % Ensure the flow is sampled at a uniform sample rate flow = retime(flow,'regular','linear','TimeStep',seconds(1e-3)); % Remove the mean from the flow and compute the flow spectrum fA = flow; fA.Data = fA.Data - mean(fA.Data); [flowSpectrum,flowFrequencies] = pspectrum(fA,'FrequencyLimits',[2 250]); % Find the values of the fault variables from the SimulationInput simin = data.SimulationInput{1}; vars = {simin.Variables.Name}; idx = strcmp(vars,'leak_cyl_area_WKSP'); LeakFault = simin.Variables(idx).Value; idx = strcmp(vars,'block_in_factor_WKSP'); BlockingFault = simin.Variables(idx).Value; idx = strcmp(vars,'bearing_fault_frict_WKSP'); BearingFault = simin.Variables(idx).Value; % Collect the fault values in a cell array faultValues = {... 'LeakFault', LeakFault, ... 'BlockingFault', BlockingFault, ... 'BearingFault', BearingFault}; end function ci = extractCI(flow,flowP,flowF) % Helper function to extract condition indicators from the flow signal % and spectrum. % Find the frequency of the peak magnitude in the power spectrum. pMax = max(flowP); fPeak = flowF(flowP==pMax); % Compute the power in the low frequency range 10-20 Hz. fRange = flowF >= 10 & flowF <= 20; pLow = sum(flowP(fRange)); % Compute the power in the mid frequency range 40-60 Hz. fRange = flowF >= 40 & flowF <= 60; pMid = sum(flowP(fRange)); % Compute the power in the high frequency range >100 Hz. fRange = flowF >= 100; pHigh = sum(flowP(fRange)); % Find the frequency of the spectral kurtosis peak [pKur,fKur] = pkurtosis(flow); pKur = fKur(pKur == max(pKur)); % Compute the flow cumulative sum range. csFlow = cumsum(flow.Data); csFlowRange = max(csFlow)-min(csFlow); % Collect the feature and feature values in a cell array. % Add flow statistic (mean, variance, etc.) and common signal % characteristics (rms, peak2peak, etc.) to the cell array. ci = {... 'qMean', mean(flow.Data), ... 'qVar', var(flow.Data), ... 'qSkewness', skewness(flow.Data), ... 'qKurtosis', kurtosis(flow.Data), ... 'qPeak2Peak', peak2peak(flow.Data), ... 'qCrest', peak2rms(flow.Data), ... 'qRMS', rms(flow.Data), ... 'qMAD', mad(flow.Data), ... 'qCSRange',csFlowRange, ... 'fPeak', fPeak, ... 'pLow', pLow, ... 'pMid', pMid, ... 'pHigh', pHigh, ... 'pKurtosis', pKur(1)}; end
参考
関連するトピック
You can also select a web site from the following list:
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
