Main Content

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

데이터 기반 모델을 사용한 결함 검출

이 예제에서는 결함 검출에 데이터 기반 모델링 접근 방식을 사용하는 방법을 보여줍니다.

소개

기계 작동의 이상 상태를 조기에 감지하고 분리하면 사고를 줄이고 가동 중단 시간을 줄여 운영 비용을 절약하는 데 도움이 됩니다. 이 접근 방식에서는 시스템 작동의 실시간 측정값을 처리하여 새로 진행된 결함을 가리킬 수 있는 예기치 않은 동작에 플래그를 지정합니다.

이 예제에서는 다음과 같은 결함 진단 측면을 살펴봅니다.

  1. 잔차 분석에 의한 이상 시스템 동작 감지

  2. 손상된 시스템의 모델을 구축하여 품질 저하 감지

  3. 모델 파라미터의 온라인 적응을 사용하여 시스템 변화 추적

시스템 동작의 동적 모델 식별하기

모델 기반 감지 접근 방식에서는 먼저 측정된 입력 및 출력 데이터를 사용하여 해당 시스템의 동적 모델을 구축합니다. 좋은 모델은 미래의 특정 시간 지평에 대해 시스템의 응답을 정확하게 예측할 수 있습니다. 예측이 좋지 않으면 잔차가 커질 수 있고 상관을 포함할 수 있습니다. 이러한 측면을 바탕으로 고장 발생 여부를 감지합니다.

충격과 진동이 가해지는 건물을 살펴보겠습니다. 진동의 원인은 돌풍, 구동 중인 엔진 및 터빈과의 접촉, 지반 진동과 같이 시스템에 따라 다양한 유형의 자극일 수 있습니다. 충격은 시스템에 더해져서 시스템에 충분한 자극을 가하는 임펄스성 범프 테스트의 결과입니다. Simulink 모델 pdmMechanicalSystem.slx는 이러한 구조의 간단한 예입니다. 가진(excitation)은 주기적 범프뿐만 아니라 필터링된 백색 잡음으로 모델링되는 지반 진동에서도 옵니다. 시스템의 출력값은 측정 잡음이 가해지는 센서에 의해 수집됩니다. 모델은 정상 상태 또는 결함 상태의 구조에 대한 다양한 시나리오를 시뮬레이션할 수 있습니다.

sysA = 'pdmMechanicalSystem';
open_system(sysA)
% Set the model in the healthy mode of operation
set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','Normal')
% Simulate the system and log the response data
sim(sysA)
ynormal = logsout.getElement('y').Values;

입력 신호는 측정되지 않았습니다. 단지 응답 ynormal만 기록되었습니다. 따라서 "블라인드 식별" 기법을 사용하여 시스템의 동적 모델을 구축합니다. 특히 여기서는 시스템을 표현하는 한 가지 방식으로서, 기록된 신호의 ARMA 모델을 구축합니다. 이 접근 방식은 입력 신호가 (필터링된) 백색 잡음이라 가정되는 경우에 유효합니다. 데이터에는 주기적 범프가 가해지므로, 데이터를 각 범프 발생 시작 시점을 기준으로 여러 조각으로 분할합니다. 이렇게 하면 각 데이터 세그먼트는 하나의 범프와 임의 가진에 대한 응답을 포함하게 됩니다. 이 상황은 시계열 모델을 사용하여 캡처할 수 있으며, 여기서 범프의 영향은 적합한 초기 상태에 기인합니다.

Ts = 1/256;  % data sample time
nr = 10;     % number of bumps in the signal
N = 512;     % length of data between bumps
znormal = cell(nr,1);
for ct = 1:nr
   ysegment = ynormal.Data((ct-1)*N+(1:500));
   z = iddata(ysegment,[],Ts);
   znormal{ct} = z;  % each segment has only one bump
end
plot(znormal{:}) % plot a sampling of the recorded segments
title('Measured Response Segments')

데이터를 추정 조각과 검증 조각으로 분할합니다.

ze = merge(znormal{1:5});
zv = merge(znormal{6:10});

ssest() 명령을 사용하여 상태공간 형식의 7차 시계열 모델을 추정합니다. 모델 차수는 교차 검증(검증 데이터에 대한 피팅 확인) 및 잔차 분석(잔차는 서로 상관관계가 없음을 확인)을 통해 선택되었습니다.

nx = 7;
model = ssest(ze, nx, 'form', 'canonical', 'Ts', Ts);
present(model)  % view model equations with parameter uncertainty
                                                                              
model =                                                                       
  Discrete-time identified state-space model:                                 
    x(t+Ts) = A x(t) + K e(t)                                                 
       y(t) = C x(t) + e(t)                                                   
                                                                              
  A =                                                                         
                       x1                  x2                  x3             
   x1                   0                   1                   0             
   x2                   0                   0                   1             
   x3                   0                   0                   0             
   x4                   0                   0                   0             
   x5                   0                   0                   0             
   x6                   0                   0                   0             
   x7  0.5548 +/- 0.04606   -2.713 +/- 0.2198    5.885 +/- 0.4495             
                                                                              
                       x4                  x5                  x6             
   x1                   0                   0                   0             
   x2                   0                   0                   0             
   x3                   1                   0                   0             
   x4                   0                   1                   0             
   x5                   0                   0                   1             
   x6                   0                   0                   0             
   x7    -8.27 +/- 0.5121    9.234 +/- 0.3513   -7.956 +/- 0.1408             
                                                                              
                       x7                                                     
   x1                   0                                                     
   x2                   0                                                     
   x3                   0                                                     
   x4                   0                                                     
   x5                   0                                                     
   x6                   1                                                     
   x7   4.263 +/- 0.02599                                                     
                                                                              
  C =                                                                         
       x1  x2  x3  x4  x5  x6  x7                                             
   y1   1   0   0   0   0   0   0                                             
                                                                              
  K =                                                                         
                      y1                                                      
   x1  1.025 +/- 0.01401                                                      
   x2  1.444 +/-  0.0131                                                      
   x3  1.907 +/- 0.01271                                                      
   x4  2.385 +/- 0.01203                                                      
   x5  2.857 +/- 0.01456                                                      
   x6   3.26 +/-  0.0222                                                      
   x7  3.552 +/-  0.0336                                                      
                                                                              
Sample time: 0.0039062 seconds                                                
                                                                              
Parameterization:                                                             
   CANONICAL form with indices: 7.                                            
   Disturbance component: estimate                                            
   Number of free coefficients: 14                                            
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.
                                                                              
Status:                                                                       
Termination condition: Near (local) minimum, (norm(g) < tol)..                
Number of iterations: 7, Number of function evaluations: 15                   
                                                                              
Estimated using SSEST on time domain data "ze".                               
Fit to estimation data: [99.07 99.04 99.15 99.05 99.04]% (prediction focus)   
FPE: 0.6242, MSE: [0.5974 0.6531 0.5991 0.5871 0.6496]                        
More information in model's "Report" property.                                
                                                                              

모델 디스플레이는 파라미터 추정값에서 상대적으로 작은 불확실성을 보입니다. 측정된 신호의 추정된 스펙트럼에 대해 1 표준편차(99.73%) 신뢰한계를 계산하여 신뢰성을 확인할 수 있습니다.

h = spectrumplot(model);
showConfidence(h, 3)

낮은 주파수에서 응답에 약 30%의 불확실성이 있긴 하지만 신뢰영역은 작습니다. 검증의 다음 단계는 모델이 검증 데이터셋 zv에서 응답을 얼마나 잘 예측하는지를 확인하는 것입니다. 25개 스텝의 사전 예측 지평을 사용합니다.

compare(zv, model, 25) % Validation against one dataset

플롯은 모델이 미래의 검증 데이터셋 25개 시간 스텝(= 0.1초)의 첫 번째 실험에서 85%가 넘는 정확도로 응답을 예측할 수 있음을 보여줍니다. 데이터셋에 있는 다른 실험에 대한 피팅을 확인하려면 플롯 좌표축을 마우스 오른쪽 버튼으로 클릭하여 나타나는 상황별 메뉴를 사용하십시오.

모델 검증의 마지막 단계는 모델에서 생성된 잔차를 분석하는 것입니다. 좋은 모델에서는 잔차가 백색입니다. 즉, 0이 아닌 지연값(lags)에 대해 통계적으로 유의미하지 않은 상관을 보입니다.

resid(model, zv)

잔차는 0이 아닌 지연값에서는 대부분 상관관계가 없습니다. 정상 동작의 모델을 도출했으니 이번에는 모델을 사용하여 결함을 검출하는 방법을 살펴봅니다.

정상 상태 모델을 사용하는 잔차 분석에 의한 결함 검출

결함 검출이란 시스템 관측값에서 원치 않거나 예기치 않은 변화에 태그를 지정하는 것입니다. 결함은 시스템 동특성의 변화를 유발하는데, 이는 점진적인 마모 및 마손에 의한 것일 수도 있고, 센서 고장이나 부서진 부품으로 인한 급격한 변화에 의한 것일 수도 있습니다. 결함이 발생할 경우, 정상 동작 상태에서 얻은 모델은 관측된 응답을 예측할 수 없습니다. 이로 인해 측정된 응답과 예측된 응답 간의 차이(잔차)가 증가합니다. 이러한 편차는 일반적으로 잔차 제곱합의 값이 크거나 상관이 존재하는 경우 플래그가 지정됩니다.

Simulink 모델을 손상된 시스템 Variant에 배치하고 시뮬레이션합니다. 잔차 테스트에는 초기 상태로 인한 과도를 보일 가능성이 있는 백색 입력값이 필요하므로 단일 범프를 입력값으로 사용합니다.

set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DamagedSystem');
set_param([sysA,'/Pulse'],'Period','5120') % to force only one bump
sim(sysA)
y = logsout.getElement('y').Values;

resid(model, y.Data)
set_param([sysA,'/Pulse'],'Period','512') % restore original

이제 잔차가 더 커졌고 0이 아닌 지연값에서 상관을 보입니다. 이처럼 잔차 메트릭을 생성하여 각각의 새로운 측정값 세트에 대해 어떻게 변화하는지를 관찰하는 것이 바로 결함 검출을 뒷받침하는 기본 아이디어입니다. 여기서는 단일 스텝 예측 오차를 기반으로 하는 간단한 잔차를 사용했습니다. 실제 현장에서는 응용 분야의 필요에 따라 맞춤 설정된 더 정교한 잔차를 생성합니다.

정상 및 품질 저하 상태의 모델을 사용한 결함 검출

결함 검출에 대한 더 세밀한 접근 방식은 시스템의 결함(손상된) 상태의 모델도 식별하는 것입니다. 그러면 어느 모델이 시스템의 라이브 측정값을 설명할 가능성이 더 큰지를 분석할 수 있습니다. 이 방식은 여러 유형의 결함에 대한 모델로 일반화될 수 있으며, 결함의 검출뿐 아니라 어느 결함인지를 식별("분리")하는 데도 사용될 수 있습니다. 이 예제에서는 다음과 같은 접근 방식을 사용합니다.

  1. 정상 상태 및 알려진 마모 및 마손으로 유도된 수명 종료 상태에서 작동하는 시스템으로부터 데이터를 수집합니다.

  2. 각 상태의 동작을 나타내는 동적 모델을 식별합니다.

  3. 데이터 군집화 접근 방식을 사용하여 여러 상태를 명확히 구분합니다.

  4. 결함 검출을 위해서, 구동 중인 기계에서 데이터를 수집하여 그 동작의 모델을 식별합니다. 그런 다음 어느 상태(정상 또는 손상)가 관측된 동작을 가장 잘 설명하는지 예측합니다.

앞에서 이미 정상 작동 모드의 시스템을 시뮬레이션했습니다. 이번에는 "수명 종료" 모드에 있는 모델 pdmMechanicalSystem을 시뮬레이션합니다. 이것은 허용 가능한 마지막 상태에 다다를 정도로 시스템이 악화된 시나리오입니다.

set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DamagedSystem');
sim(sysA)
y = logsout.getElement('y').Values;
zfault = cell(nr,1);
for ct = 1:nr
   z = iddata(y.Data((ct-1)*N+(1:500)),[],Ts);
   zfault{ct} = z;
end

각 데이터 세그먼트에 대해 모델 세트를 만듭니다. 앞에서와 같이 상태공간 형식의 7차 시계열 모델을 구축합니다. 속도를 높이기 위해 공분산 계산은 해제합니다.

mNormal =  cell(nr,1);
mFault = cell(nr, 1);
nx = order(model);
opt = ssestOptions('EstimateCovariance',0);
for ct = 1:nr
   mNormal{ct} = ssest(znormal{ct}, nx, 'form', 'canonical', 'Ts', Ts, opt);
   mFault{ct} = ssest(zfault{ct}, nx, 'form', 'canonical', 'Ts', Ts, opt);
end

모델 mFault가 결함 작동 모드를 효과적으로 표현하는지 다음처럼 확인합니다.

compare(merge(zfault{:}), mFault{:}, 25)

아래에 정상 추정 스펙트럼과 결함 추정 스펙트럼이 플로팅되어 있습니다.

Color1 = 'k'; Color2 = 'r';
ModelSet1 = cat(2,mNormal,repmat({Color1},[nr, 1]))';
ModelSet2 = cat(2,mFault,repmat({Color2},[nr, 1]))';

spectrum(ModelSet1{:},ModelSet2{:})
axis([1  1000  -45  40])
title('Output Spectra (black: normal, red: faulty)')

스펙트럼 플롯은 두 모델 간의 차이를 보여줍니다. 손상된 모드의 주 공명이 증폭되지만 그 밖의 경우에는 두 스펙트럼이 중첩됩니다. 다음으로, 정상 상태와 결함 상태를 정량적으로 구분할 방법을 만듭니다. 다음과 같은 데이터 군집화 및 분류 접근 방식을 사용할 수 있습니다.

  • 퍼지 C-평균 군집화 Fuzzy Logic Toolbox의 fcm()을 참조하십시오.

  • 서포트 벡터 머신 분류기. Statistics and Machine Learning Toolbox의 fitcsvm ()을 참조하십시오.

  • 자기 조직화 맵. Deep Learning Toolbox의 selforgmap()을 참조하십시오.

이 예제에서는 서포트 벡터 머신 분류 기법을 사용합니다. 두 유형의 모델(mNormalmFault)에서 얻어진 정보의 군집화는 이들 모델이 제공할 수 있는 여러 유형의 정보(예: 모델의 극점과 영점의 위치, 모델의 피크 공명의 위치, 모델의 파라미터 목록)를 기반으로 할 수 있습니다. 여기서는 두 공명에 대응되는 극점 위치를 기준으로 모드를 분류합니다. 군집화의 경우, 정상 상태 모델의 극점에 'good'이라는 태그를 지정하고 결함 상태 모델의 극점에 'faulty'라는 태그를 지정합니다.

ModelTags = cell(nr*2,1);  % nr is number of data segments
ModelTags(1:nr) = {'good'};
ModelTags(nr+1:end) = {'faulty'};
ParData = zeros(nr*2,4);
plist = @(p)[real(p(1)),imag(p(1)),real(p(3)),imag(p(3))]; % poles of dominant resonances
for ct = 1:nr
   ParData(ct,:) =  plist(esort(pole(mNormal{ct})));
   ParData(nr+ct,:) = plist(esort(pole(mFault{ct})));
end
cl = fitcsvm(ParData,ModelTags,'KernelFunction','rbf', ...
   'BoxConstraint',Inf,'ClassNames',{'good', 'faulty'});
cl.ConvergenceInfo.Converged
ans =

  logical

   1

cl은 훈련 데이터 ParData를 정상 영역과 결함 영역으로 분리하는 SVM 분류기입니다. 이 분류기의 predict 메서드를 사용하면 입력 nx×1 벡터를 두 영역 중 하나에 할당할 수 있습니다.

이제 분류기의 예측(정상 대 손상)을 테스트하고, 시스템이 연속적으로 정상(모드 = 'Normal')에서 완전 손상(모드 = 'DamagedSystem')으로 변해가는 방식으로 파라미터가 변화하는 시스템으로부터 데이터 배치를 수집할 수 있습니다. 이 시나리오를 시뮬레이션하기 위해 모델을 'DeterioratingSystem' 모드에 둡니다.

set_param([sysA,'/Mechanical System'],'LabelModeActiveChoice','DeterioratingSystem');
sim(sysA)
ytv = logsout.getElement('y').Values; ytv = squeeze(ytv.Data);
PredictedMode = cell(nr,1);
ssOpt = ssestOptions('SearchMethod','lm');
for ct = 1:nr
   zSegment = iddata(ytv((ct-1)*512+(1:500)),[],Ts);
   mSegment = ssest(zSegment, nx, 'form', 'canonical', 'Ts', Ts, ssOpt);
   PredictedMode(ct) = predict(cl, plist(esort(pole(mSegment))));
end

I = strcmp(PredictedMode,'good');
Tags = ones(nr,1);
Tags(~I) = -1;
t = (0:5120)'*Ts;  % simulation time
Time = t(1:512:end-1);
plot(Time(I),Tags(I),'g*',Time(~I),Tags(~I),'r*','MarkerSize',12)
grid on
axis([0 20 -2 2])
title('Green: Normal, Red: Faulty state')
xlabel('Data evaluation time')
ylabel('Prediction')

플롯은 분류기가 대략 중간 점까지는 동작을 정상으로 예측하고 이후부터는 결함 상태로 예측함을 보여줍니다.

모델 파라미터의 온라인 적응에 의한 결함 검출

앞에 나온 분석에서는 시스템 작동 중 여러 시점에서 수집한 데이터 배치를 사용했습니다. 시스템의 건전성을 모니터링하는 보다 편리한 대안으로는 시스템 동작의 적응 모델을 만드는 것이 있습니다. 새로운 측정값은 지속적으로 처리되며 모델의 파라미터를 재귀적으로 업데이트하는 데 사용됩니다. 마모 및 마손 또는 결함의 영향은 모델 파라미터 값의 변화로 나타납니다.

마모 및 마손 시나리오를 다시 살펴보겠습니다. 시스템이 노후화됨에 따라 "덜컹거림"이 심해지고, 이는 여러 공명 모드의 가진과 시스템 피크 응답의 증가로 나타납니다. 이 시나리오는 모델 pdmDeterioratingSystemEstimation에 설명되어 있습니다. 이 모델은 오프라인 식별을 위해 추가된 임펄스성 범프가 존재하지 않는다는 점을 제외하면 pdmMechanicalSystem의 'DeterioratingSystem' 모드와 동일합니다. 시스템의 응답은 ARMA 모델 구조의 파라미터를 추정하도록 구성된 "Recursive Polynomial Model Estimator" 블록으로 전달됩니다. 실제 시스템은 정상 상태에서 시작하지만, 200초의 시간 범위에 걸쳐 수명 종료 상태로 품질이 저하됩니다.

initial_model = translatecov(@(x)idpoly(x),model);
sysB = 'pdmDeterioratingSystemEstimation';
open_system(sysB);

"ARMA 모델" 블록은 이전 섹션에서 도출한 정상 동작의 추정 모델에서 얻어진 파라미터 및 공분산 데이터를 다항식(ARMA) 형식으로 변환한 후에 사용하여 초기화되었습니다. 파라미터 공분산 데이터도 변환될 수 있도록 translatecov() 함수가 사용되었습니다. 이 블록은 "망각 인자(Forgetting factor)" 알고리즘을 사용하며, 망각 인자는 각 샘플링 시점에 파라미터를 업데이트하도록 1보다 약간 작은 값으로 설정되었습니다. 망각 인자의 값은 시스템이 얼마나 빠르게 업데이트되는지에 영향을 줍니다. 값이 작으면 업데이트의 분산이 커지고, 값이 크면 추정기가 빠른 변화에 적응하기가 어려워집니다.

모델 파라미터 추정값은 출력 스펙트럼과 출력 스펙트럼의 3 표준편차 신뢰영역을 업데이트하는 데 사용됩니다. 스펙트럼의 신뢰영역이 관심 주파수에서 정상 시스템의 신뢰영역과 중첩되지 않으면 시스템이 변한 것을 명확히 알 수 있습니다. 플롯에 검은색 선으로 표시된 결함 검출 임계값은 특정 주파수에서 허용되는 최대 이득을 나타냅니다. 시스템의 변화가 누적됨에 따라 스펙트럼이 이 선을 넘나듭니다. 이는 결함에 대한 시각적 지표로 쓰일 수 있으며, 실제 시스템에서는 이를 사용하여 수리를 요청할 수 있습니다.

시뮬레이션을 실행하고 스펙트럼 플롯이 업데이트되는 모습을 살펴봅니다.

sim(sysB)

모델 파라미터의 실행 추정값은 시스템 극점 위치를 계산하는 데에도 사용되며, 극점 위치는 SVM 분류기에 입력되어 시스템이 "정상" 상태와 "결함" 상태 중 어느 상태에 있는지 예측하는 데 사용됩니다. 이 결정도 플롯에 표시되어 있습니다. 예측의 정규화된 점수가 0.3보다 작으면 결정은 잠정적인 것으로(구분 경계에 가까움) 간주됩니다. 스펙트럼 및 분류기 예측의 실행 추정값이 계산되는 방식은 스크립트 pdmARMASpectrumPlot.m을 참조하십시오.

Simulink 외부에서 recursiveARMA() 함수를 사용하여 적응 추정 및 플로팅 절차를 구현하는 것이 가능합니다. "Recursive Polynomial Model Estimator" 블록과 recursiveARMA() 함수 모두 배포를 위한 코드 생성을 지원합니다.

분류 방식은 알려진 고장 모드가 여러 개인 경우로 일반화될 수 있습니다. 이를 위해서는 하나의 모드가 특정한 유형의 고장을 나타내는 다중 그룹 분류기가 필요합니다. 이 예제에서는 이러한 경우를 다루지 않습니다.

결론

이 예제에서는 시스템 식별 방식에 데이터 군집화 및 분류 접근 방식을 결합하면 결함의 검출 및 분리에 어떤 도움이 되는지 보여주었습니다. 순차 배치 분석 방식과 온라인 적응 방식을 모두 살펴보았습니다. 측정된 출력 신호의 ARMA 구조 모델이 식별되었습니다. 입력 신호와 출력 신호를 모두 사용할 수 있고 상태공간 모델, Box-Jenkins 다항식 모델과 같은 다른 유형의 모델 구조를 사용하려는 경우에도 비슷한 접근 방식을 사용할 수 있습니다.

이 예제에서는 다음 내용을 확인했습니다.

  1. 정상 작동 모델을 기반으로 하는 잔차의 상관은 고장이 시작되었음을 나타낼 수 있습니다.

  2. 점진적으로 악화되는 결함은 시스템 동작의 지속적인 적응 모델을 사용하여 검출할 수 있습니다. 모델 특성에 대한 사전 설정 임계값(예: 출력 스펙트럼의 범위)은 고장의 시작 및 진행을 시각화하는 데 도움이 될 수 있습니다.

  3. 결함의 원인을 분리해야 하는 경우, 사전에 해당 고장 모드에 대해 별도의 모델을 만드는 접근 방식을 사용할 수 있습니다. 그런 다음 분류 접근 방식을 사용하여 시스템의 예측된 상태를 모델 중 하나에 할당할 수 있습니다.

관련 항목