잔차 분석을 사용한 원심 펌프의 결함 진단

다음 제품이 필요합니다.

이 예제에서는 모델 패리티 방정식 기반의 방식을 사용하여 펌핑 시스템에서 발생하는 여러 유형의 결함을 검출하고 진단하는 방법을 보여줍니다. 이 예제에서는 정상 상태 실험을 사용한 원심 펌프의 결함 진단 예제에서 제시한 기법을 펌프 작동이 복수의 작동 상태에서 이루어지는 상황으로 확장해서 적용해 봅니다.

이 예제는 Rolf Isermann의 저서 Fault Diagnosis Applications에서 제시하는 원심 펌프 분석을 따릅니다 [1]. System Identification Toolbox™, Statistics and Machine Learning Toolbox™, Control System Toolbox™, Simulink™의 기능을 사용합니다. Predictive Maintenance Toolbox™는 필요하지 않습니다.

다중 속도 펌프 구동 - 잔차 분석에 의한 진단

정상 상태 펌프의 헤드 및 토크 방정식은 펌프가 빠르게 변하는 속도로 구동되거나 넓은 범위의 속도에서 구동되는 경우에는 정확한 결과를 생성하지 않습니다. 마찰을 비롯한 기타 손실이 유의미해질 수 있으며, 모델의 파라미터는 속도에 대한 종속성을 보입니다. 이러한 경우에 널리 적용할 수 있는 접근 방식은 동작의 블랙박스 모델을 만드는 것입니다. 이러한 모델의 파라미터는 물리적으로 유의미할 필요가 없습니다. 모델은 알려진 동작의 시뮬레이션을 위한 장치로 사용됩니다. 모델의 출력값을 대응되는 측정된 신호에서 뺄셈하여 잔차를 계산합니다. 잔차의 속성(평균, 분산, 전력 등)은 정상 작동과 결함 작동을 구분하는 데 사용됩니다.

정적 펌프 헤드 방정식과 동적 펌프-파이프 방정식을 사용하여, 그림에 나와 있는 4개의 잔차를 계산할 수 있습니다.

모델에는 다음과 같은 컴포넌트가 있습니다.

정적 펌프 모델: $Δ p_{}^{ˆ} (t) = θ_{1} ω^{2} (t) + θ_{2} ω (t)$
동적 파이프 모델: $Q_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{Δ p (t)} + θ_{5} Q_{}^{ˆ} (t - 1)$
동적 펌프-파이프 모델: ${Q_{}^{ˆ}}_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{{Δ p (t)}_{}^{ˆ}} + θ_{5} {Q_{}^{ˆ}}_{}^{ˆ} (t - 1)$
동적 역 펌프 모델: ${M_{}^{ˆ}}_{m o t o r} (t) = θ_{6} + θ_{7} ω (t) + θ_{8} ω (t - 1) + θ_{9} ω^{2} (t) + θ_{10} {M_{}^{ˆ}}_{m o t o r} (t - 1)$

모델 파라미터 $θ_{1}, . . ., θ_{10}$ 은 펌프 속도에 대한 종속성을 보입니다. 이 예제에서는 파라미터의 조각별 선형 근사가 계산됩니다. 작동 영역을 3개의 영역으로 나눕니다.

$ω \leq 900 R P M$
$900 < ω \leq 1500 R P M$
$ω > 1500 R P M$

정상 펌프를 폐루프 제어기를 갖는 폐루프에서 기준 속도 범위 0~3000RPM으로 구동했습니다. 기준 입력은 수정된 PRBS 신호입니다. 모터 토크, 펌프 토크, 속도 및 압력의 측정값은 10Hz 샘플링 주파수에서 수집되었습니다. 측정된 신호를 불러오고 기준 펌프 속도와 실제 펌프 속도를 플로팅합니다. 이 데이터는 최대 20MB의 대용량 데이터셋으로 MathWorks 지원 파일 사이트에서 다운로드됩니다.

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/DynamicOperationData.mat';
websave('DynamicOperationData.mat',url);
load DynamicOperationData
figure
plot(t, RefSpeed, t, w)
xlabel('Time (s)')
ylabel('Pump Speed (RPM)')
legend('Reference','Actual')

펌프 속도 범위를 기반으로 작동 영역을 정의합니다.

I1 = w<=900;            % first operating regime
I2 = w>900 & w<=1500;   % second operating regime
I3 = w>1500;            % third operating regime

모델 식별

A. 정적 펌프 모델 식별

펌프 속도 $ω (t)$ 와 차동 압력 $Δ p (t)$ 의 측정된 값을 입력-출력 데이터로 사용하여 정적 펌프 방정식의 파라미터 $θ_{1}$ 과 $θ_{2}$ 를 추정합니다. 이 추정을 수행하는 헬퍼 함수 staticPumpEst를 참조하십시오.

th1 = zeros(3,1);  
th2 = zeros(3,1); 
dpest = nan(size(dp));  % estimated pressure difference
[th1(1), th2(1), dpest(I1)] = staticPumpEst(w, dp, I1);  % Theta1, Theta2 estimates for regime 1
[th1(2), th2(2), dpest(I2)] = staticPumpEst(w, dp, I2);  % Theta1, Theta2 estimates for regime 2
[th1(3), th2(3), dpest(I3)] = staticPumpEst(w, dp, I3);  % Theta1, Theta2 estimates for regime 3
plot(t, dp, t, dpest) % compare measured and predicted pressure differential
xlabel('Time (s)')
ylabel('\Delta P')
legend('Measured','Estimated','Location','best')
title('Static Pump Model Validation')

B. 동적 파이프 모델 식별

유량 $Q (t)$ 와 차동 압력 $Δ p (t)$ 의 측정된 값을 입력-출력 데이터로 사용하여 파이프 배출 흐름 방정식 $Q_{}^{ˆ} (t) = θ_{3} + θ_{4} \sqrt{Δ p (t)} + θ_{5} Q_{}^{ˆ} (t - 1)$ 의 파라미터 $θ_{3}$ , $θ_{4}$ , $θ_{5}$ 를 추정합니다. 이 추정을 수행하는 헬퍼 함수 dynamicPipeEst를 참조하십시오.

th3 = zeros(3,1);  
th4 = zeros(3,1); 
th5 = zeros(3,1);
[th3(1), th4(1), th5(1)] = dynamicPipeEst(dp, Q, I1); % Theta3, Theta4, Theta5 estimates for regime 1
[th3(2), th4(2), th5(2)] = dynamicPipeEst(dp, Q, I2); % Theta3, Theta4, Theta5 estimates for regime 2
[th3(3), th4(3), th5(3)] = dynamicPipeEst(dp, Q, I3); % Theta3, Theta4, Theta5 estimates for regime 3

동적 파이프 모델은 정적 펌프 모델의 경우와 달리 유량 값에 대한 동적 종속성을 보입니다. 가변 속도 영역에서 모델을 시뮬레이션하기 위해 Simulink에서 Control System Toolbox의 "LPV System" 블록을 사용하여 조각별 선형 모델을 생성합니다. Simulink 모델 LPV_pump_pipe와 이 시뮬레이션을 수행하는 헬퍼 함수 simulatePumpPipeModel을 참조하십시오.

% Check Control System Toolbox availability
ControlsToolboxAvailable = ~isempty(ver('control')) && license('test', 'Control_Toolbox');
if ControlsToolboxAvailable
    % Simulate the dynamic pipe model. Use measured value of pressure as input
    Ts = t(2)-t(1);
    Switch = ones(size(w));
    Switch(I2) = 2;
    Switch(I3) = 3;
    UseEstimatedP = 0;
    Qest_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);
    plot(t,Q,t,Qest_pipe) % compare measured and predicted flow rates
else
    % Load pre-saved simulation results from the piecewise linear Simulink model
    load DynamicOperationData Qest_pipe
    Ts = t(2)-t(1);
    plot(t,Q,t,Qest_pipe)
end
xlabel('Time (s)')
ylabel('Flow rate (Q), m^3/s')
legend('Measured','Estimated','Location','best')
title('Dynamic Pipe Model Validation')

C. 동적 펌프 파이프 모델 식별

동적 펌프-파이프 모델은 위에서 식별된 동일한 파라미터( $θ_{3}, θ_{4}, θ_{5}$ )를 사용하되, 모델 시뮬레이션에 측정된 압력 차이가 아닌 추정된 압력 차이가 필요하다는 점만 다릅니다. 따라서 새로운 식별이 필요하지 않습니다. $θ_{3}, θ_{4}, θ_{5}$ 의 추정값이 펌프-파이프 동특성을 양호하게 재현하는지 확인합니다.

if ControlsToolboxAvailable
    UseEstimatedP = 1;
    Qest_pump_pipe = simulatePumpPipeModel(Ts,th3,th4,th5);
    plot(t,Q,t,Qest_pump_pipe) % compare measured and predicted flow rates
else
    load DynamicOperationData Qest_pump_pipe 
    plot(t,Q,t,Qest_pump_pipe)
end

xlabel('Time (s)')
ylabel('Flow rate Q (m^3/s)')
legend('Measured','Estimated','location','best')
title('Dynamic Pump-Pipe Model Validation')

피팅은 측정된 압력 값을 사용하여 얻은 것과 실질적으로 동일합니다.

D. 동적 역 펌프 모델 식별

파라미터 $θ_{6}, . . ., θ_{10}$ 은 이전의 토크 및 속도 측정값에 대해 측정된 토크 값을 회귀하여 비슷한 방법으로 식별할 수 있습니다. 그러나 결과로 생성되는 조각별 선형 모델의 완전한 다중 속도 시뮬레이션은 데이터에 대한 양호한 피팅을 제공하지 않습니다. 따라서 유리 회귀 변수를 사용하여 데이터를 피팅하는 비선형 ARX 모델의 식별이 수반되는 또 다른 블랙박스 모델링 접근 방식이 사용됩니다.

% Use first 300 samples out of 550 for identification
N = 350;
sys3 = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)

sys3 = 
Nonlinear ARX model with 1 output and 2 inputs
  Inputs: u1, u2
  Outputs: y1

Regressors:
  1. Linear regressors in variables y1, u1, u2
  2. Custom regressor: u1(t-2)^2
  3. Custom regressor: u1(t)*u2(t-2)
  4. Custom regressor: u2(t)^2
  List of all regressors

Output function: Linear Function
Sample time: 0.1 seconds

Status:                                         
Estimated using NLARX on time domain data.      
Fit to estimation data: 49.2% (simulation focus)
FPE: 1798, MSE: 3392

Mmot_est = sim(sys3,[w Q]);
plot(t,Mmot,t,Mmot_est) % compare measured and predicted motor torque
xlabel('Time (s)')
ylabel('Motor Torque (Nm)')
legend('Measured','Estimated','location','best')
title('Inverse pump model validation')

잔차 생성

모델의 잔차를 측정된 신호와 모델에서 생성된 이에 대응되는 출력값 사이의 차이로 정의합니다. 네 개의 모델 컴포넌트에 대응되는 네 개의 잔차를 계산합니다.

r1 = dp - dpest;
r2 = Q - Qest_pipe;
r3 = Q - Qest_pump_pipe;

역 펌프 모델 잔차의 계산을 위해, 원래 잔차는 분산이 크므로 이동평균 필터를 사용하여 모델 출력값에 평활화 연산을 적용합니다.

r4 = Mmot - movmean(Mmot_est,[1 5]);

다음과 같이 훈련 잔차를 확인합니다.

figure
subplot(221)
plot(t,r1)
ylabel('Static pump - r1')
subplot(222)
plot(t,r2)
ylabel('Dynamic pipe - r2')
subplot(223)
plot(t,r3)
ylabel('Dynamic pump-pipe - r3')
xlabel('Time (s)')
subplot(224)
plot(t,r4)
ylabel('Dynamic inverse pump - r4')
xlabel('Time (s)')

잔차 특징 추출

잔차는 결함 분리를 위해 적합한 특징이 추출되는 신호입니다. 파라미터 정보가 없으므로, 신호의 최대 진폭이나 분산과 같은 신호 속성에서만 도출된 특징을 고려합니다.

PRBS 입력 실현을 사용하여 펌프-파이프 시스템에 수행한 20개의 실험 세트를 살펴봅니다. 실험 세트는 다음의 각 모드에 대해 반복됩니다.

정상 펌프
결함 1: 조립 틈새의 마모
결함 2: 임펠러 출구의 작은 침전
결함 3: 임펠러 주입구의 침전
결함 4: 임펠러 출구의 연삭 마모
결함 5: 부러진 블레이드
결함 6: 공동 현상
결함 7: 속도 센서 편향
결함 8: 유량계 편향
결함 9: 압력 센서 편향

실험 데이터를 불러옵니다.

url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/MultiSpeedOperationData.mat';
websave('MultiSpeedOperationData.mat',url);
load MultiSpeedOperationData
% Generate operation mode labels
Labels = {'Healthy','ClearanceGapWear','ImpellerOutletDeposit',...
    'ImpellerInletDeposit','AbrasiveWear','BrokenBlade','Cavitation','SpeedSensorBias',...
    'FlowmeterBias','PressureSensorBias'};

각 앙상블과 각 작동 모드에 대해 잔차를 계산합니다. 이 작업은 몇 분 정도 걸릴 수 있습니다. 잔차 데이터가 데이터 파일에 저장됩니다. helperComputeEnsembleResidues를 실행하여 잔차를 생성합니다.

% HealthyR = helperComputeEnsembleResidues(HealthyEnsemble,Ts,sys3,th1,th2,th3,th4,th5); % Healthy data residuals

% Load pre-saved data from "helperComputeEnsembleResidues" run
url = 'https://www.mathworks.com/supportfiles/predmaint/fault-diagnosis-of-centrifugal-pumps-using-residual-analysis/Residuals.mat';
websave('Residuals.mat',url);
load Residuals

모드 변별력이 가장 큰 잔차의 특징은 사전에 알려져 있지 않습니다. 따라서 각 잔차에 대해 평균, 최대 진폭, 분산, 첨도, 1-노름과 같은 여러 후보 특징을 생성합니다. 모든 특징은 정상 앙상블의 값 범위를 사용하여 스케일링됩니다.

CandidateFeatures = {@mean, @(x)max(abs(x)), @kurtosis, @var, @(x)sum(abs(x))};
FeatureNames = {'Mean','Max','Kurtosis','Variance','OneNorm'};
% generate feature table from gathered residuals of each fault mode
[HealthyFeature, MinMax] = helperGenerateFeatureTable(HealthyR, CandidateFeatures, FeatureNames);
Fault1Feature  = helperGenerateFeatureTable(Fault1R,  CandidateFeatures, FeatureNames, MinMax);
Fault2Feature  = helperGenerateFeatureTable(Fault2R,  CandidateFeatures, FeatureNames, MinMax);
Fault3Feature  = helperGenerateFeatureTable(Fault3R,  CandidateFeatures, FeatureNames, MinMax);
Fault4Feature  = helperGenerateFeatureTable(Fault4R,  CandidateFeatures, FeatureNames, MinMax);
Fault5Feature  = helperGenerateFeatureTable(Fault5R,  CandidateFeatures, FeatureNames, MinMax);
Fault6Feature  = helperGenerateFeatureTable(Fault6R,  CandidateFeatures, FeatureNames, MinMax);
Fault7Feature  = helperGenerateFeatureTable(Fault7R,  CandidateFeatures, FeatureNames, MinMax);
Fault8Feature  = helperGenerateFeatureTable(Fault8R,  CandidateFeatures, FeatureNames, MinMax);
Fault9Feature  = helperGenerateFeatureTable(Fault9R,  CandidateFeatures, FeatureNames, MinMax);

각 특징 테이블에는 20개의 특징이 있습니다(각 잔차 신호에 대해 5개의 특징). 각 테이블은 실험 하나당 하나씩 총 50개의 관측값(행)을 포함합니다.

N = 50; % number of experiments in each mode
FeatureTable = [...
   [HealthyFeature(1:N,:), repmat(Labels(1),[N,1])];...
   [Fault1Feature(1:N,:),  repmat(Labels(2),[N,1])];...
   [Fault2Feature(1:N,:),  repmat(Labels(3),[N,1])];...
   [Fault3Feature(1:N,:),  repmat(Labels(4),[N,1])];...
   [Fault4Feature(1:N,:),  repmat(Labels(5),[N,1])];...
   [Fault5Feature(1:N,:),  repmat(Labels(6),[N,1])];...
   [Fault6Feature(1:N,:),  repmat(Labels(7),[N,1])];...
   [Fault7Feature(1:N,:),  repmat(Labels(8),[N,1])];...
   [Fault8Feature(1:N,:),  repmat(Labels(9),[N,1])];...
   [Fault9Feature(1:N,:),  repmat(Labels(10),[N,1])]];
FeatureTable.Properties.VariableNames{end} = 'Condition';

% Preview some samples of training data
disp(FeatureTable([2 13 37 49 61 62 73 85 102 120],:))

     Mean1       Mean2      Mean3      Mean4      Max1        Max2        Max3        Max4       Kurtosis1    Kurtosis2    Kurtosis3    Kurtosis4     Variance1    Variance2    Variance3    Variance4    OneNorm1    OneNorm2    OneNorm3    OneNorm4            Condition        
    ________    _______    _______    _______    _______    ________    ________    _________    _________    _________    _________    __________    _________    _________    _________    _________    ________    ________    ________    ________    _________________________

     0.32049     0.6358     0.6358    0.74287    0.39187     0.53219     0.53219     0.041795      0.18957    0.089492     0.089489              0     0.45647       0.6263       0.6263      0.15642      0.56001    0.63541     0.63541      0.24258    {'Healthy'              }
     0.21975    0.19422    0.19422    0.83704          0     0.12761     0.12761      0.27644      0.18243     0.19644      0.19643        0.20856    0.033063      0.16772      0.16773     0.025119     0.057182    0.19344     0.19344     0.063167    {'Healthy'              }
      0.1847    0.25136    0.25136    0.87975    0.32545     0.13929      0.1393     0.084162      0.72346    0.091488     0.091485       0.052552    0.063757      0.18371      0.18371     0.023845      0.10671    0.25065     0.25065      0.04848    {'Healthy'              }
           1          1          1          0    0.78284     0.94535     0.94535       0.9287      0.41371     0.45927      0.45926        0.75934     0.92332            1            1      0.80689      0.99783          1           1       0.9574    {'Healthy'              }
     -2.6545    0.90555    0.90555    0.86324     1.3037      0.7492      0.7492      0.97823    -0.055035     0.57019      0.57018         1.4412      1.4411      0.73128      0.73128      0.80573       3.2084    0.90542     0.90542      0.49257    {'ClearanceGapWear'     }
    -0.89435    0.43384    0.43385     1.0744    0.82197     0.30254     0.30254    -0.022325      0.21411    0.098504       0.0985      -0.010587     0.55959      0.29554      0.29554     0.066693       1.0432    0.43326     0.43326      0.13785    {'ClearanceGapWear'     }
     -1.2149    0.53579    0.53579     1.0899    0.87439     0.47339     0.47339      0.29547     0.058946    0.048138     0.048137        0.17864     0.79648      0.48658      0.48658      0.25686       1.4066     0.5353      0.5353      0.26663    {'ClearanceGapWear'     }
     -1.0949    0.44616    0.44616      1.143    0.56693     0.19719     0.19719    -0.012055     -0.10819     0.32603      0.32603     -0.0033592     0.43341      0.12857      0.12857     0.038392       1.2238    0.44574     0.44574     0.093243    {'ClearanceGapWear'     }
     -0.1844    0.16651    0.16651    0.87747    0.25597    0.080265    0.080265      0.49715       0.5019     0.29939      0.29938         0.5338    0.072397     0.058024     0.058025       0.1453      0.20263    0.16566     0.16566      0.14216    {'ImpellerOutletDeposit'}
     -3.0312     0.9786     0.9786    0.75241      1.473     0.97839     0.97839       1.0343    0.0050647      0.4917       0.4917        0.89759      1.7529       0.9568       0.9568       0.8869       3.6536    0.97859     0.97859      0.62706    {'ImpellerOutletDeposit'}

분류기 설계

A. 산점도 플롯을 사용하여 모드 분리도 시각화하기

특징을 시각적으로 조사하는 것으로 분석을 시작합니다. 이를 위해 결함 1: 조립 틈새의 마모를 살펴봅니다. 이 결함을 검출하는 데 어느 특징이 가장 적합한지 확인하기 위해 'Healthy' 레이블과 'ClearanceGapWear' 레이블을 사용하여 특징의 산점도 플롯을 생성합니다.

T = FeatureTable(:,1:20);
P = T.Variables;
R = FeatureTable.Condition;
I = strcmp(R,'Healthy') | strcmp(R,'ClearanceGapWear');
f = figure;
gplotmatrix(P(I,:),[],R(I))
f.Position(3:4) = f.Position(3:4)*1.5;

뚜렷하게 나타나지는 않지만, 1열과 17열의 특징이 가장 큰 분리도를 보입니다. 이들 특징을 더 상세히 분석합니다.

f = figure;
Names = FeatureTable.Properties.VariableNames;
J = [1 17];
fprintf('Selected features for clearance gap fault: %s\n',strjoin(Names(J),', '))

Selected features for clearance gap fault: Mean1, OneNorm1

gplotmatrix(P(I,[1 17]),[],R(I))

이제 플롯은 정상 모드와 조립 틈새 결함 모드를 분리하는 데 특징 Mean1과 OneNorm1을 사용할 수 있음을 분명히 보여줍니다. 각 결함 모드에 대해서도 유사한 분석을 수행할 수 있습니다. 모든 경우에, 결함 모드를 구분하는 특징 세트를 찾을 수 있습니다. 따라서 결함 동작의 검출은 항상 가능합니다. 그러나 동일한 특징이 여러 결함 유형의 영향을 받으므로 결함 분리는 이보다 더 어렵습니다. 예를 들어, 특징 Mean1(r1의 평균)과 OneNorm1(r1의 1-노름)은 여러 결함 유형에 대해 변화를 보입니다. 그래도 센서 편향과 같은 일부 결함은 여러 특징에서 분리 가능하므로 비교적 쉽게 분리할 수 있습니다.

세 개의 센서 편향 결함에 대해, 산점도 플롯을 직접 조사하여 특징을 선택합니다.

figure;
I = strcmp(R,'Healthy') | strcmp(R,'PressureSensorBias') | strcmp(R,'SpeedSensorBias') | strcmp(R,'FlowmeterBias');
J = [1 4 6 16 20]; % selected feature indices
fprintf('Selected features for sensors'' bias: %s\n',strjoin(Names(J),', '))

Selected features for sensors' bias: Mean1, Mean4, Max2, Variance4, OneNorm4

gplotmatrix(P(I,J),[],R(I))

선택된 특징의 산점도 플롯은 하나 이상의 특징 쌍에서 4개 모드를 구분할 수 있음을 보여줍니다. 축소된 특징 세트를 사용하여 축소된 결함 세트(센서 편향 결함만)에 대해 20개 멤버의 TreeBagger 분류기를 훈련시킵니다.

rng default % for reproducibility
Mdl = TreeBagger(20, FeatureTable(I,[J 21]), 'Condition',...
   'OOBPrediction','on',...
   'OOBPredictorImportance','on');
figure
plot(oobError(Mdl))
xlabel('Number of trees')
ylabel('Misclassification probability')

오분류 오차가 3% 미만입니다. 따라서 특정 결함 서브셋의 분류를 위해 작은 특징 세트를 선택하여 사용하는 것이 가능합니다.

B. 분류 학습기 앱을 사용한 다중 클래스 분류

이전 섹션에서는 산점도 플롯을 직접 조사하여 특정 결함 유형에 대해 특징 세트를 축소해 보았습니다. 이 접근 방식은 번거로울 수 있으며 모든 결함 유형을 다루기 어려울 수 있습니다. 보다 자동화된 방식으로 모든 결함 모드를 처리하는 분류기를 설계할 수 있을까요? Statistics and Machine Learning Toolbox에는 다양한 분류기가 있습니다. 분류 학습기 앱을 사용하면 여러 분류기를 빠르게 사용해 보고 성능을 비교해 볼 수 있습니다.

분류 학습기 앱을 실행하고, 새로운 세션에서 사용할 작업 데이터로 작업 공간에 있는 FeatureTable을 선택합니다. 데이터의 20%(각 모드당 10개 샘플)는 홀드아웃 검증을 위해 남겨 둡니다.
메인 탭의 모델 유형 섹션에서 모두를 선택합니다. 그런 다음 훈련 버튼을 누릅니다.
짧은 시간 안에 20개의 분류기가 훈련됩니다. 내역 패널에서 분류기 이름 옆에 정확도가 표시됩니다. 선형 SVM 분류기가 홀드아웃 샘플에 대해 86%의 정확도를 보이므로 가장 성능이 좋습니다. 이 분류기는 "ClearanceGapWear"를 분류하는 데 얼마간의 어려움을 보이는데, 이 결함을 40%의 확률로 "ImpellerOutletDeposit"으로 분류합니다.
성능을 그래픽으로 보려면 메인 탭의 플롯 섹션에서 혼동행렬 앱을 여십시오. 플롯은 선택한 분류기(여기서는 선형 SVM 분류기)의 성능을 보여줍니다.

가장 성능이 좋은 분류기를 작업 공간으로 내보내고 이를 새로운 측정값에 대한 예측에 사용합니다.

요약

결함 진단 전략을 잘 설계하면 서비스 가동 중단 시간 및 부품 교체 비용을 최소화하여 운영 비용을 절감할 수 있습니다. 작동하는 기계의 동특성에 대한 지식과 센서 측정값을 함께 사용하여 여러 종류의 결함을 검출하고 분리하면 훌륭한 전략을 수립할 수 있습니다. 이 예제에서는 원심 펌프의 결함 진단을 위한 잔차 기반 접근 방식을 설명했습니다. 이 접근 방식은 모델링 작업이 복잡하고 모델 파라미터가 작동 상태에 대한 종속성을 보이는 경우 파라미터 추정 및 추적 기반 접근 방식의 좋은 대안이 됩니다.

잔차 기반 결함 진단 접근 방식에는 다음과 같은 단계가 사용됩니다.

물리적 고려 사항 또는 블랙박스 시스템 식별 기법을 사용하여 시스템의 측정 가능한 입력값과 출력값 사이의 동특성을 모델링합니다.
잔차를 측정된 신호와 모델에서 생성된 신호 사이의 차이로 계산합니다. 결함 분리도를 개선하려면 잔차를 추가로 필터링해야 할 수 있습니다.
각 잔차 신호에서 피크 진폭, 전력, 첨도 등과 같은 특징을 추출합니다.
이상 감지 및 분류 기법을 사용하여 결함 검출 및 분류를 위해 특징을 사용합니다.
잔차와 도출된 특징 전부가 모든 결함에 민감한 것은 아닙니다. 특징 히스토그램 및 산점도 플롯 보기에서 특정 결함 유형을 검출하는 데 어느 특징이 적합한지 확인할 수 있습니다. 특징을 선택하고 결함 분리 성능을 평가하는 과정은 반복적인 절차일 수 있습니다.

참고 문헌

Isermann, Rolf, Fault-Diagnosis Applications. Model-Based Condition Monitoring: Actuators, Drives, Machinery, Plants, Sensors, and Fault-tolerant System, Edition 1, Springer-Verlag Berlin Heidelberg, 2011.

지원 함수

정적 펌프 방정식 파라미터 추정

function [x1, x2, dpest] = staticPumpEst(w, dp, I)
%staticPumpEst Static pump parameter estimation in a varying speed setting
% I: sample indices for the selected operating region.

w1 = [0; w(I)];
dp1 = [0; dp(I)];
R1 = [w1.^2 w1];
x = pinv(R1)*dp1;
x1 = x(1);  
x2 = x(2);  

dpest = R1(2:end,:)*x;
end

동적 파이프 파라미터 추정

function [x3, x4, x5, Qest] = dynamicPipeEst(dp, Q, I)
%dynamicPipeEst Dynamic pipe parameter estimation in a varying speed setting
% I: sample indices for the selected operating region.

Q = Q(I);
dp = dp(I);
R1 = [0; Q(1:end-1)];
R2 = dp; R2(R2<0) = 0; R2 = sqrt(R2);
R = [ones(size(R2)), R2, R1];

% Remove out-of-regime samples
ii = find(I);
j = find(diff(ii)~=1);
R = R(2:end,:); R(j,:) = [];
y = Q(2:end); y(j) = [];
x = R\y;

x3 = x(1);
x4 = x(2);
x5 = x(3);

Qest = R*x;
end

LPV System 블록을 사용하여 펌프-파이프 모델의 동적 다중 작동 모드를 시뮬레이션합니다.

function Qest = simulatePumpPipeModel(Ts,th3,th4,th5)
%simulatePumpPipeModel Piecewise linear modeling of dynamic pipe system.
% Ts: sample time
% w: Pump rotational speed
% th1, th2, th3 are estimated model parameters for the 3 regimes.
% This function requires Control System Toolbox.

ss1 = ss(th5(1),th4(1),th5(1),th4(1),Ts);
ss2 = ss(th5(2),th4(2),th5(2),th4(2),Ts);
ss3 = ss(th5(3),th4(3),th5(3),th4(3),Ts);
offset = permute([th3(1),th3(2),th3(3)]',[3 2 1]);
OP = struct('Region',[1 2 3]');
sys = cat(3,ss1,ss2,ss3);
sys.SamplingGrid = OP;

assignin('base','sys',sys)
assignin('base','offset',offset)
mdl = 'LPV_pump_pipe';
sim(mdl);
Qest = logsout.get('Qest');
Qest = Qest.Values;
Qest = Qest.Data;
end

역 펌프 동특성의 동적 모델을 식별합니다.

function syse = identifyNonlinearARXModel(Mmot,w,Q,Ts,N)
%identifyNonlinearARXModel Identify a nonlinear ARX model for 2-input (w, Q), 1-output (Mmot) data.
% Inputs:
%  w: rotational speed
%  Q: Flow rate
%  Mmot: motor torque
%  N: number of data samples to use
% Outputs:
%  syse: Identified model
%
% This function uses NLARX estimator from System Identification Toolbox.

sys = idnlarx([2 2 1 0 1],'','CustomRegressors',{'u1(t-2)^2','u1(t)*u2(t-2)','u2(t)^2'});
data = iddata(Mmot,[w Q],Ts);
opt = nlarxOptions;
opt.Focus = 'simulation';
opt.SearchOptions.MaxIterations = 500;
syse = nlarx(data(1:N),sys,opt);
end