확장 칼만 필터를 사용한 결함 검출

라이브 스크립트 열기

이 예제에서는 결함 검출을 위해 확장 칼만 필터를 사용하는 방법을 보여줍니다. 이 예제에서는 간단한 DC 모터의 마찰을 온라인 추정하기 위해 확장 칼만 필터를 사용합니다. 추정된 마찰에 유의미한 변화가 감지되고 이는 결함을 나타냅니다. 이 예제에서는 System Identification Toolbox™의 기능을 사용하며, Predictive Maintenance Toolbox™는 필요하지 않습니다.

모터 모델

모터는 토크 u로 구동되며 감쇠 계수 c를 갖는 관성 J로 모델링됩니다. 모터 각속도 w와 가속도 $w_{}^{˙}$ 는 측정된 출력값입니다.

$w_{}^{˙} = (u - c w) / J$

확장 칼만 필터를 사용하여 감쇠 계수 c를 추정하기 위해 감쇠 계수의 보조 상태를 도입하고 그 도함수를 0으로 설정합니다.

$c_{}^{˙} = 0$

따라서 모델 상태 x = [w;c]와 측정값 y의 방정식은 다음과 같습니다.

$[\begin{array}{c} w_{}^{˙} \\ c_{}^{˙} \end{array}] = [\begin{array}{c} (u - c w) / J \\ 0 \end{array}]$

$y = [\begin{array}{c} w \\ (u - c w) / J \end{array}]$

연속시간 방정식은 근사 $x_{}^{˙} = \frac{x_{n + 1} - x_{n}}{T_{s}}$ 을 사용하여 이산시간으로 변환됩니다. 여기서 Ts는 이산 샘플링 주기입니다. 이 이산시간 모델 방정식은 pdmMotorModelStateFcn.m 함수와 pdmMotorModelMeasurementFcn.m 함수에 구현됩니다.

$[\begin{array}{c} w_{n + 1} \\ c_{n + 1} \end{array}] = [\begin{array}{c} w_{n} + (u_{n} - c_{n} w_{n}) T_{s} / J \\ c_{n} \end{array}]$

$y_{n} = [\begin{array}{c} w_{n} \\ (u_{n} - c_{n} w_{n}) / J \end{array}]$

모터 파라미터를 지정합니다.

J  = 10;    % Inertia 
Ts = 0.01;  % Sample time

초기 상태를 지정합니다.

x0 = [...
    0; ...  % Angular velocity 
    1];     % Friction

type pdmMotorModelStateFcn

function x1 = pdmMotorModelStateFcn(x0,varargin)
%PDMMOTORMODELSTATEFCN
%
% State update equations for a motor with friction as a state
%
%  x1 = pdmMotorModelStateFcn(x0,u,J,Ts)
%
%  Inputs:
%    x0 - initial state with elements [angular velocity; friction] 
%    u  - motor torque input
%    J  - motor inertia
%    Ts - sampling time
%
%  Outputs:
%    x1 - updated states
%

%  Copyright 2016-2017 The MathWorks, Inc.

% Extract data from inputs
u  = varargin{1};   % Input
J  = varargin{2};   % System innertia
Ts = varargin{3};   % Sample time

% State update equation
x1 = [...
    x0(1)+Ts/J*(u-x0(1)*x0(2)); ...
    x0(2)];
end

type pdmMotorModelMeasurementFcn

function y = pdmMotorModelMeasurementFcn(x,varargin)
%PDMMOTORMODELMEASUREMENTFCN
%
% Measurement equations for a motor with friction as a state
%
%  y = pdmMotorModelMeasurementFcn(x0,u,J,Ts)
%
%  Inputs:
%    x  - motor state with elements [angular velocity; friction] 
%    u  - motor torque input
%    J  - motor inertia
%    Ts - sampling time
%
%  Outputs:
%    y - motor measurements with elements [angular velocity; angular acceleration]
%

%  Copyright 2016-2017 The MathWorks, Inc.

% Extract data from inputs
u  = varargin{1};   % Input
J  = varargin{2};   % System innertia

% Output equation
y = [...
    x(1); ...
    (u-x(1)*x(2))/J];
end

모터에는 상태(공정) 잡음 외란 q와 측정 잡음 외란 r이 가해집니다. 잡음 항은 가산성입니다.

$[\begin{array}{c} w_{n + 1} \\ c_{n + 1} \end{array}] = [\begin{array}{c} w_{n} + (u_{n} - c_{n} w_{n}) T_{s} / J \\ c_{n} \end{array}] + q$

$y_{n} = [\begin{array}{c} w_{n} \\ (u_{n} - c_{n} w_{n}) / J \end{array}] + r$

공정 및 측정 잡음은 평균 0 E[q]=E[r]=0과 공분산 Q = E[qq'] 및 R = E[rr']을 갖습니다. 마찰 상태는 높은 공정 잡음 외란을 갖습니다. 이는 모터가 정상적으로 작동하는 동안 마찰에 변동이 있을 것으로 예상하며 필터가 이 변동을 추적하기를 원한다는 점을 반영합니다. 가속도 및 속도 상태 잡음은 낮지만 속도 및 가속도 측정값은 상대적으로 잡음이 많습니다.

공정 잡음 공분산을 지정합니다.

Q = [...
    1e-6 0; ...   % Angular velocity
    0 1e-2];      % Friction

측정 잡음 공분산을 지정합니다.

R = [...
    1e-4 0; ...  % Velocity measurement
    0 1e-4];     % Acceleration measurement

확장 칼만 필터 생성하기

모델의 상태를 추정하기 위해 확장 칼만 필터를 만듭니다. 감쇠 상태 값의 급격한 변화는 결함 이벤트를 나타내므로 감쇠 상태에 특히 관심이 있습니다.

extendedKalmanFilter 객체를 만들고 상태 천이 및 측정 함수의 야코비 행렬을 지정합니다.

ekf = extendedKalmanFilter(...
    @pdmMotorModelStateFcn, ...
    @pdmMotorModelMeasurementFcn, ...
    x0,...
    'StateCovariance',            [1 0; 0 1000], ...[1 0 0; 0 1 0; 0 0 100], ...
    'ProcessNoise',               Q, ...
    'MeasurementNoise',           R, ...
    'StateTransitionJacobianFcn', @pdmMotorModelStateJacobianFcn, ...
    'MeasurementJacobianFcn',     @pdmMotorModelMeasJacobianFcn);

확장 칼만 필터는 앞에서 정의된 상태 천이 및 측정 함수를 입력 인수로 갖습니다. 초기 상태 값 x0, 초기 상태 공분산, 공정 및 측정 잡음 공분산도 확장 칼만 필터의 입력값입니다. 이 예제에서는 상태 천이 함수 f와 측정 함수 h에서 정확한 야코비 행렬 함수를 도출할 수 있습니다.

$\frac{\partial}{\partial x} f = [\begin{array}{cc} 1 - T_{s} c_{n} / J & - T_{s} w_{n} / J \\ 0 & 1 \end{array}]$

$\frac{\partial}{\partial x} h = [\begin{array}{cc} 1 & 0 \\ - c_{n} / J & - w_{n} / J \end{array}]$

상태 야코비 행렬은 pdmMotorModelStateJacobianFcn.m 함수에 정의되어 있고 측정 야코비 행렬은 pdmMotorModelMeasJacobianFcn.m 함수에 정의되어 있습니다.

type pdmMotorModelStateJacobianFcn

function Jac = pdmMotorModelStateJacobianFcn(x,varargin)
%PDMMOTORMODELSTATEJACOBIANFCN
%
% Jacobian of motor model state equations. See pdmMotorModelStateFcn for
% the model equations.
%
%  Jac = pdmMotorModelJacobianFcn(x,u,J,Ts)
%
%  Inputs:
%    x  - state with elements [angular velocity; friction] 
%    u  - motor torque input
%    J  - motor inertia
%    Ts - sampling time
%
%  Outputs:
%    Jac - state Jacobian computed at x
%

%  Copyright 2016-2017 The MathWorks, Inc.

% Model properties
J  = varargin{2};
Ts = varargin{3};

% Jacobian
Jac = [...
    1-Ts/J*x(2) -Ts/J*x(1); ...
    0 1];
end

type pdmMotorModelMeasJacobianFcn

function J = pdmMotorModelMeasJacobianFcn(x,varargin)
%PDMMOTORMODELMEASJACOBIANFCN
%
% Jacobian of motor model measurement equations. See
% pdmMotorModelMeasurementFcn for the model equations.
%
%  Jac = pdmMotorModelMeasJacobianFcn(x,u,J,Ts)
%
%  Inputs:
%    x  - state with elements [angular velocity; friction] 
%    u  - motor torque input
%    J  - motor inertia
%    Ts - sampling time
%
%  Outputs:
%    Jac - measurement Jacobian computed at x
%

%  Copyright 2016-2017 The MathWorks, Inc.

% System parameters
J  = varargin{2};   % System innertia

% Jacobian
J = [ ...
    1 0;
    -x(2)/J -x(1)/J];
end

시뮬레이션

플랜트를 시뮬레이션하기 위해 루프를 만들고 모터에 결함(모터 마찰의 급격한 변화)을 발생시킵니다. 시뮬레이션 루프 안에서 확장 칼만 필터를 사용하여 모터 상태를 추정하고, 특히 마찰 상태를 추적하여 마찰에 통계적으로 유의미한 변화가 발생하는 시점을 감지합니다.

모터는 모터를 반복적으로 가속 및 감속하는 펄스 열을 사용하여 시뮬레이션됩니다. 이러한 유형의 모터 작동은 생산 라인의 선택기 로봇에서 전형적으로 볼 수 있습니다.

t  = 0:Ts:20;                  % Time, 20s with Ts sampling period
u  = double(mod(t,1)<0.5)-0.5; % Pulse train, period 1, 50% duty cycle
nt = numel(t);                 % Number of time points
nx = size(x0,1);               % Number of states
ySig = zeros([2, nt]);         % Measured motor outputs
xSigTrue = zeros([nx, nt]);    % Unmeasured motor states
xSigEst = zeros([nx, nt]);     % Estimated motor states
xstd = zeros([nx nx nt]);      % Standard deviation of the estimated states 
ySigEst = zeros([2, nt]);      % Estimated model outputs
fMean = zeros(1,nt);           % Mean estimated friction
fSTD = zeros(1,nt);            % Standard deviation of estimated friction 
fKur = zeros(2,nt);            % Kurtosis of estimated friction
fChanged = false(1,nt);        % Flag indicating friction change detection

모터를 시뮬레이션할 때는 확장 칼만 필터를 생성할 때 사용한 Q 및 R 잡음 공분산 값과 비슷한 공정 및 측정 잡음을 추가합니다. 마찰의 경우, 마찰은 결함이 발생하는 경우를 제외하고 대부분 일정하므로 훨씬 작은 잡음 값을 사용합니다. 시뮬레이션 중에 인공적으로 결함을 유발합니다.

rng('default');
Qv = chol(Q);   % Standard deviation for process noise
Qv(end) = 1e-2; % Smaller friction noise
Rv = chol(R);   % Standard deviation for measurement noise

상태 업데이트 방정식을 사용하여 모델을 시뮬레이션하고, 모델 상태에 공정 잡음을 추가합니다. 시뮬레이션을 시작한 지 10초가 지나면 모터 마찰에 강제로 변화를 줍니다. 모델 측정 함수를 사용하여 모터 센서를 시뮬레이션하고, 모델 출력값에 측정 잡음을 추가합니다.

for ct = 1:numel(t)
    
   % Model output update   
   y = pdmMotorModelMeasurementFcn(x0,u(ct),J,Ts);
   y = y+Rv*randn(2,1);   % Add measurement noise
   ySig(:,ct) = y;
   
   % Model state update 
   xSigTrue(:,ct) = x0;
   x1 = pdmMotorModelStateFcn(x0,u(ct),J,Ts);
   % Induce change in friction
   if t(ct) == 10
       x1(2) = 10;  % Step change
   end
   x1n = x1+Qv*randn(nx,1);  % Add process noise
   x1n(2) = max(x1n(2),0.1); % Lower limit on friction
   x0 = x1n; % Store state for next simulation iteration

모터 측정값으로부터 모터 상태를 추정하기 위해 확장 칼만 필터의 predict 명령과 correct 명령을 사용합니다.

   % State estimation using the Extended Kalman Filter
   x_corr = correct(ekf,y,u(ct),J,Ts); % Correct the state estimate based on current measurement.
   xSigEst(:,ct) = x_corr;
   xstd(:,:,ct) = chol(ekf.StateCovariance);
   predict(ekf,u(ct),J,Ts);            % Predict next state given the current state and input.

마찰의 변화를 감지하기 위해, 추정된 마찰의 평균 및 표준편차를 4초 이동 윈도우를 사용하여 계산합니다. 첫 7초가 지난 후에, 계산된 평균 및 표준편차를 고정합니다. 처음에 계산된 이 평균은 마찰에 대해 예상되는 결함 없음(no-fault) 평균값입니다. 7초가 지난 후에, 추정된 마찰이 예상 결함 없음 평균값에서 3 표준편차보다 더 많이 떨어져 있다면 이는 마찰의 유의미한 변화를 나타냅니다. 추정된 마찰에 포함된 잡음과 변동성의 영향을 줄이고자, 3 표준편차 범위와 비교할 때는 추정된 마찰의 평균을 사용합니다.

   if t(ct) < 7 
       % Compute mean and standard deviation of estimated fiction.
       idx = max(1,ct-400):max(1,ct-1); % Ts = 0.01 seconds
       fMean(ct) = mean( xSigEst(2, idx) );
       fSTD(ct)  = std( xSigEst(2, idx) );
   else
       % Store the computed mean and standard deviation without
       % recomputing.
       fMean(ct) = fMean(ct-1);
       fSTD(ct)  = fSTD(ct-1);
       % Use the expected friction mean and standard deviation to detect
       % friction changes.
       estFriction = mean(xSigEst(2,max(1,ct-10):ct));
       fChanged(ct) = (estFriction > fMean(ct)+3*fSTD(ct)) || (estFriction < fMean(ct)-3*fSTD(ct));
   end
   if fChanged(ct) && ~fChanged(ct-1) 
       % Detect a rising edge in the friction change signal |fChanged|.
       fprintf('Significant friction change at %f\n',t(ct));
   end

Significant friction change at 10.450000

추정된 상태를 사용하여 추정된 출력값을 계산합니다. 측정된 출력값과 추정된 출력값 사이의 오차를 계산하고, 오차 통계량을 계산합니다. 오차 통계량은 마찰 변화를 감지하는 데 사용할 수 있습니다. 이 부분은 뒤에서 더 자세히 설명합니다.

   ySigEst(:,ct) = pdmMotorModelMeasurementFcn(x_corr,u(ct),J,Ts);   
   idx = max(1,ct-400):ct;
   fKur(:,ct) = [...
       kurtosis(ySigEst(1,idx)-ySig(1,idx)); ...
       kurtosis(ySigEst(2,idx)-ySig(2,idx))];
end

확장 칼만 필터의 성능

10.45초에서 마찰 변화가 감지되었다는 점에 유의하십시오. 이번에는 이 결함 검출 규칙이 어떻게 도출되었는지 설명합니다. 먼저 시뮬레이션 결과와 필터 성능을 검토합니다.

figure, 
subplot(211), plot(t,ySig(1,:),t,ySig(2,:));
title('Motor Outputs')
legend('Measured Angular Velocity','Measured Angular Acceleration', 'Location','SouthWest')
subplot(212), plot(t,u);
title('Motor Input - Torque')

Figure contains 2 axes objects. Axes object 1 with title Motor Outputs contains 2 objects of type line. These objects represent Measured Angular Velocity, Measured Angular Acceleration. Axes object 2 with title Motor Input - Torque contains an object of type line.

모델의 입력-출력 응답을 보면 측정된 신호에서 직접 마찰 변화를 감지하기 어려운 것을 알 수 있습니다. 확장 칼만 필터를 사용하면 상태, 그중에서도 특히 마찰 상태를 추정할 수 있습니다. 실제 모델 상태와 추정된 상태를 비교합니다. 추정된 상태를 3 표준편차에 해당하는 신뢰구간과 함께 표시했습니다.

figure, 
subplot(211),plot(t,xSigTrue(1,:), t,xSigEst(1,:), ...
    [t nan t],[xSigEst(1,:)+3*squeeze(xstd(1,1,:))', nan, xSigEst(1,:)-3*squeeze(xstd(1,1,:))'])
axis([0 20 -0.06 0.06]), 
legend('True value','Estimated value','Confidence interval')
title('Motor State - Velocity')
subplot(212),plot(t,xSigTrue(2,:), t,xSigEst(2,:),  ...
    [t nan t],[xSigEst(2,:)+3*squeeze(xstd(2,2,:))' nan xSigEst(2,:)-3*squeeze(xstd(2,2,:))'])
axis([0 20 -10 15])
title('Motor State - Friction');

Figure contains 2 axes objects. Axes object 1 with title Motor State - Velocity contains 3 objects of type line. These objects represent True value, Estimated value, Confidence interval. Axes object 2 with title Motor State - Friction contains 3 objects of type line.

필터로 추정된 값은 실제 값을 추종하고, 신뢰구간은 유계로 유지되는 것을 볼 수 있습니다. 추정 오차를 검토하면 필터 동작을 더 상세히 살펴볼 수 있습니다.

figure, 
subplot(211),plot(t,xSigTrue(1,:)-xSigEst(1,:))
title('Velocity State Error')
subplot(212),plot(t,xSigTrue(2,:)-xSigEst(2,:))
title('Friction State Error')

Figure contains 2 axes objects. Axes object 1 with title Velocity State Error contains an object of type line. Axes object 2 with title Friction State Error contains an object of type line.

오차 플롯은 10초에서 발생하는 마찰 변화 후에 필터가 적응하여 추정 오차를 0으로 줄임을 보여줍니다. 그러나 오차 플롯은 실제 상태를 알고 있어야 도출할 수 있으므로 결함 검출에 사용할 수 없습니다. 가속도 및 속도의 측정된 상태 값을 추정된 상태 값과 비교하면 검출 메커니즘을 도출하는 것이 가능할 수 있습니다.

figure
subplot(211), plot(t,ySig(1,:)-ySigEst(1,:))
title('Velocity Measurement Error')
subplot(212),plot(t,ySig(2,:)-ySigEst(2,:))
title('Acceleration Measurement Error')

Figure contains 2 axes objects. Axes object 1 with title Velocity Measurement Error contains an object of type line. Axes object 2 with title Acceleration Measurement Error contains an object of type line.

가속도 오차 플롯은 결함이 발생한 10초 부근에서 평균 오차에 약간 차이가 있음을 보여줍니다. 오차 통계량을 검토하여, 계산된 오차에서 결함을 검출할 수 있는지 살펴봅니다. 가속도 및 속도 오차는 정규분포를 따를 것으로 예상됩니다(잡음 모델은 모두 가우스 분포입니다). 따라서 가속도 오차의 첨도는 오차 분포가 마찰 변화 및 그로 인한 오차 분포의 변화로 인해 대칭에서 비대칭으로 바뀌는 시점을 식별하는 데 도움이 될 수 있습니다.

figure,
subplot(211),plot(t,fKur(1,:))
title('Velocity Error Kurtosis')
subplot(212),plot(t,fKur(2,:))
title('Acceleration Error Kurtosis')

Figure contains 2 axes objects. Axes object 1 with title Velocity Error Kurtosis contains an object of type line. Axes object 2 with title Acceleration Error Kurtosis contains an object of type line.

추정기가 아직 수렴 중이고 데이터가 수집되고 있는 처음 4초를 무시하면, 오차의 첨도는 3(가우스 분포의 예상 첨도 값) 부근에서 약간 변동이 있으면서 비교적 일정한 것을 볼 수 있습니다. 따라서 오차 통계량은 이 응용 분야에서 마찰 변화를 자동으로 감지하는 데 사용할 수 없습니다. 이 응용 분야에서는 오차의 첨도를 사용하기도 어렵습니다. 필터가 적응하고 있고 지속적으로 오차를 0에 가깝게 만들고 있어서 오차 분포가 0이 아닌 시간 윈도우가 매우 짧기 때문입니다.

따라서 이 응용 분야에서는 추정된 마찰의 변화를 사용하는 것이 모터의 결함을 자동으로 검출하는 가장 좋은 방법이 됩니다. 알려진 결함 없음 데이터에서 얻어진 마찰 추정값(평균 및 표준편차)은 마찰에 대한 예상 범위를 제공하며, 범위를 위반한 경우 쉽게 감지할 수 있습니다. 다음 플롯에서는 이 결함 검출 접근 방식을 보여줍니다.

figure
plot(t,xSigEst(2,:),[t nan t],[fMean+3*fSTD,nan,fMean-3*fSTD])
title('Friction Change Detection')
legend('Estimated Friction','No-Fault Friction Bounds')
axis([0 20 -10 20])
grid on

Figure contains an axes object. The axes object with title Friction Change Detection contains 2 objects of type line. These objects represent Estimated Friction, No-Fault Friction Bounds.

요약

이 예제에서는 확장 칼만 필터를 사용하여 간단한 DC 모터의 결함을 추정하고 결함 검출을 위해 마찰 추정값을 사용하는 방법을 보여주었습니다.

확장 칼만 필터를 사용한 결함 검출

모터 모델

확장 칼만 필터 생성하기

시뮬레이션

확장 칼만 필터의 성능

요약

관련 항목