Main Content

이 번역 페이지는 최신 내용을 담고 있지 않습니다. 최신 내용을 영문으로 보려면 여기를 클릭하십시오.

시뮬레이션된 데이터를 사용한 다중 클래스 결함 검출

이 예제에서는 Simulink 모델을 사용하여 결함 데이터와 정상 데이터를 생성하는 방법을 보여줍니다. 이 데이터를 사용하여 여러 결함 조합을 검출하는 다중 클래스 분류기를 개발합니다. 이 예제에서는 삼중 왕복 펌프 모델을 사용하며, 누수, 막힘, 베어링 결함을 포함합니다.

모델 설정

이 예제에서는 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'])

펌프 모델은 실린더 누수, 주입구 막힘, 베어링 마찰 증가와 같은 세 가지 결함 유형을 모델링하도록 구성됩니다. 이들 결함은 작업 공간 변수로 파라미터화되고 펌프 블록 대화 상자를 통해 구성됩니다.

결함 데이터와 정상 데이터 시뮬레이션하기

세 가지 결함 유형 각각에 대해 결함 없음부터 상당한 결함까지의 범위를 갖는, 결함 심각도를 나타내는 값으로 구성된 배열을 만듭니다.

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

펌프 모델은 잡음을 포함하도록 구성되므로 동일한 결함 파라미터 값으로 모델을 실행해도 서로 다른 시뮬레이션 출력값이 생성됩니다. 이는 동일한 결함 상태 및 심각도에 대해 여러 개의 시뮬레이션 결과가 생성될 수 있음을 의미하므로 분류기를 개발할 때 유용합니다. 이러한 결과를 위한 시뮬레이션을 구성하려면 값이 결함 없음, 단일 결함, 두 가지 결함의 조합, 세 가지 결함의 조합을 나타내는 결함 파라미터 값으로 구성된 벡터를 만드십시오. 각 그룹(결함 없음, 단일 결함 등)에 대해 이렇게 정의한 결함 파라미터 값으로부터 125개의 결함 값 조합을 만듭니다. 이를 통해 총 1,000개의 결함 파라미터 값의 조합이 만들어집니다. 1,000개의 시뮬레이션을 병렬로 실행하면 표준 데스크탑에서 약 1시간이 소요되며 약 620MB의 데이터가 생성됩니다. 시뮬레이션 시간을 줄이려면 runAll = truerunAll = false. Note that a로 변경하여 오류 조합의 개수를 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')

유속 스펙트럼은 예상한 주파수에 공명 피크가 있음을 나타냅니다. 구체적으로 보면, 펌프 모터 속도는 950rpm/15.833Hz이고, 펌프에 실린더가 3개 있으므로 유속은 3*15.833Hz=47.5Hz에 기본주파수를, 47.5Hz의 배수에 고조파를 가질 것으로 예상됩니다. 유속 스펙트럼에서 예상되는 공명 피크를 명확히 확인할 수 있습니다. 펌프의 실린더 중 하나에서 발생한 결함은 펌프 모터 속도 15.833Hz와 그 고조파에서 공명으로 나타납니다.

유속 스펙트럼과 느린 신호를 통해 상태 지표가 될 수 있을 만한 특징을 대략적으로 파악할 수 있습니다. 구체적으로는 평균, 분산과 같은 일반적인 신호 통계량과 스펙트럼 특징이 여기에 해당합니다. 피크 크기를 갖는 주파수, 15.833Hz 부근의 에너지, 47.5Hz 부근의 에너지, 100Hz 이상에서의 에너지와 같은 예상되는 고조파와 관련된 스펙트럼 상태 지표가 계산됩니다. 스펙트럼 첨도 피크의 주파수도 계산됩니다.

상태 지표를 위한 데이터 변수와 결함 변수 값을 위한 상태 변수로 앙상블을 구성합니다. 그런 다음 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 명령을 사용하여 앙상블의 모든 멤버에 대해 이 작업을 반복합니다. 여러 결함 상태의 시뮬레이션 결과를 파악하기 위해 앙상블의 백 번째 요소마다 플로팅합니다.

%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 

각 앙상블 멤버의 결함 변수 값(데이터 테이블에서 행으로 나타남)을 여러 결함 플래그로 변환할 수 있으며, 여러 결함 플래그는 각 멤버의 서로 다른 결함 상태를 캡처하는 하나의 플래그로 결합할 수 있습니다.

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

정오 플롯은 각 결함 조합에 대해 결함 조합이 올바르게 예측된 횟수(플롯에서 대각선에 있는 항목)와 결함 조합이 잘못 예측된 횟수(대각선을 벗어난 항목)를 보여줍니다.

정오 플롯은 분류기가 일부 결함 상태를 올바르게 분류하지 않았음을 보여줍니다(대각선을 벗어난 항목). 그러나 결함 없음 상태는 올바르게 예측되었습니다. 두 군데에서 결함이 있음에도 결함 없음 상태가 예측되었지만(첫 번째 열), 그 외의 경우에는 정확한 결함 상태는 아니더라도 결함이 예측되기는 했습니다. 전체적으로 검증 정확도는 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 

참고 항목

관련 항목