シミュレーションデータを使用した複数クラス故障検出

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

この例では次を使用します:

この例では、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 の組み合わせが得られます。

nPerGroup = 125; % Number of elements in each fault group
    
% 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 を作成します。

これら 1000 のシミュレーションを並列で実行する場合、標準のデスクトップでは約 1 時間かかり、約 620 MB のデータが生成されることに注意してください。便宜上、最初の 10 個のシミュレーションのみを実行するオプションが用意されています。

% Run the simulation and create an ensemble to manage the simulation results
runAll = true;
if runAll
    [ok,e] = generateSimulationEnsemble(simInput,fullfile('.','Data'),'UseParallel',true);
else
    [ok,e] = generateSimulationEnsemble(simInput(1:10),fullfile('.','Data')); %#ok<UNRCH>
end

[09-Apr-2018 09:01:38] Checking for availability of parallel pool...
[09-Apr-2018 09:01:38] Loading Simulink on parallel workers...
Analyzing and transferring files to the workers ...done.
[09-Apr-2018 09:01:38] Configuring simulation cache folder on parallel workers...
[09-Apr-2018 09:01:38] Running SetupFcn on parallel workers...
[09-Apr-2018 09:01:39] Loading model on parallel workers...
[09-Apr-2018 09:01:39] Transferring base workspace variables used in the model to parallel workers...
[09-Apr-2018 09:01:41] Running simulations...
[09-Apr-2018 09:02:28] Completed 1 of 1000 simulation runs
[09-Apr-2018 09:02:33] Completed 2 of 1000 simulation runs
[09-Apr-2018 09:02:37] Completed 3 of 1000 simulation runs
[09-Apr-2018 09:02:41] Completed 4 of 1000 simulation runs
[09-Apr-2018 09:02:46] Completed 5 of 1000 simulation runs
[09-Apr-2018 09:02:49] Completed 6 of 1000 simulation runs
[09-Apr-2018 09:02:54] Completed 7 of 1000 simulation runs
[09-Apr-2018 09:02:58] Completed 8 of 1000 simulation runs
[09-Apr-2018 09:03:01] Completed 9 of 1000 simulation runs...

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

シミュレーション結果からの特徴の処理と抽出

このモデルはポンプの吐出圧力、出流量、モーター速度、およびモーター電流をログに記録するよう設定されます。

ens.DataVariables

ans = 8×1 string array
    "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));

Starting parallel pool (parpool) using the 'local' profile ...
Preserving jobs with IDs: 1 2 because they contain crash dump files.
You can use 'delete(myCluster.Jobs)' to remove all jobs created with profile local. To create 'myCluster' use 'myCluster = parcluster('local')'.
connected to 6 workers.
Evaluating tall expression using the Parallel Pool 'local':
- Pass 1 of 1: Completed in 45 sec
Evaluation completed in 45 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.86472    117.63    18.874     276.49      35.572    7.5242    -0.72832      2.7738        13.835      1.1494    35.677    2.2326     42690         1e-09          0.8              0      
    43.909    0.44477    125.92    18.899     12.417      35.576     7.869     -0.7094      2.6338        13.335      1.1449    35.686    2.3204     42697         1e-09          0.8              0      
    43.909     1.1782    137.99    17.526     11.589      35.573    7.4367    -0.72208      2.7136        12.641      1.1395    35.678    2.2407     42695         1e-09          0.8              0      
    44.151     156.74    173.84    21.073      199.5      33.768    12.466    -0.30256      2.4782        17.446      1.2138    33.952    2.8582     40518       2.4e-06          0.8              0      
    43.848    0.71756    110.92    18.579     197.02      35.563    7.5781    -0.72377       2.793         14.14      1.1504    35.669    2.2671     42682         1e-09          0.8              0      
    43.909    0.43673    119.56    20.003     11.589       35.57    7.5028    -0.74797      2.7913        13.833      1.1551    35.676    2.2442     42689         1e-09          0.8              0      
    43.788    0.31617     135.3    19.724     476.82      35.568    7.4406    -0.70964      2.6884        14.685      1.1473    35.673    2.2392     42687         1e-09          0.8              0      
    43.848    0.72747    121.63    19.733     11.589      35.523     7.791    -0.72736      2.7864        14.043      1.1469    35.633    2.2722     42633         1e-09          0.8              0      
    43.848    0.62777    128.85    19.244     11.589      35.541    7.5698     -0.6953      2.6942        13.451      1.1415    35.647    2.2603     42654         1e-09          0.8              0      
    43.848     0.4631    134.83    18.918     12.417      35.561    7.8607    -0.68417      2.6664        13.885      1.1504    35.671    2.3078     42681         1e-09          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.6150

% Compute the accuracy of the classifier at predicting 
% that there is a fault
1-sum(confdata(2:end,1))/sum(confdata(:))

ans = 0.9450

故障なしと予測され、故障が存在したケースについて調べます。まず、検証データで、実際の故障が閉塞故障でありながら、故障なしと予測されたケースを見つけます。

vData = data(cvp.test(1),:);
b1 = (actualValue==2) & (predictedValue==0);
fData = vData(b1,15:17)

fData=11×3 table
    LeakFault    BlockingFault    BearingFault
    _________    _____________    ____________

      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.71              0      
      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.62              0      
      1e-09          0.77              0      
      1e-09          0.77              0      
      1e-09          0.71              0      
      8e-07          0.74              0      
      1e-09          0.74              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: Completed in 38 sec
Evaluation completed in 38 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

参考

simulationEnsembleDatastore