拡張カルマンフィルターを使用した故障検出

ライブスクリプトを開く

この例は、故障検出のために拡張カルマンフィルターを使用する方法を示します。この例では、拡張カルマンフィルターを使用して簡単な DC モーターの摩擦のオンライン推定を行います。推定された摩擦に大きな変化が検出されます。これは故障を示します。この例では、System Identification Toolbox™ の機能を使用し、Predictive Maintenance Toolbox™ は必要ありません。

モーターモデル

モーターはトルク u で駆動される、減衰係数 c をもつ慣性 J としてモデル化されます。モーターの角速度 w および角加速度 $w_{}^{˙}$ は測定された出力です。

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

拡張カルマンフィルターを使って減衰係数 c を推定するには、減衰係数に補助状態を導入してその微係数をゼロに設定します。

$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$

プロセスノイズと測定ノイズは、ゼロ平均 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 秒が経過したら、計算された平均と標準偏差をロックします。この初期計算で求めた平均は、非故障状態にある摩擦の期待される平均値です。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 秒の時点における摩擦の変化以後適応し、推定誤差をゼロに減らしていることがわかります。しかし、誤差のプロットは真の状態に関する情報に依存するため、故障検出に使用することはできません。加速度および速度について、測定された状態値と推定された状態値を比較することで、検出メカニズムが得られる可能性があります。

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 (ガウス分布で予想される尖度値) で比較的一定しています。したがって、このアプリケーションでは摩擦の変化を誤り統計によって自動的に検出することはできません。このアプリケーションでは、フィルターが適応して誤差を常にゼロにし、誤差分布がゼロ以外の値になる時間枠が短くなるため、誤差の尖度を使用するのは困難です。

そのため、このアプリケーションでは、推定された摩擦の変化を使用することがモーターの故障を自動的に検出する最良の方法です。既知の非故障時データから得た摩擦の推定値 (平均と標準偏差) に基づいて摩擦の期待範囲がわかるので、これらの範囲を超えた場合は簡単に検出できます。次のプロットは、この故障検出手法をよく示しています。

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 モーターの摩擦を推定し、推定された摩擦を使って故障を検出する方法を説明しました。

拡張カルマン フィルターを使用した故障検出

モーター モデル

拡張カルマン フィルターの作成

シミュレーション

拡張カルマン フィルターの性能

まとめ

関連するトピック

拡張カルマンフィルターを使用した故障検出

モーターモデル

拡張カルマンフィルターの作成

拡張カルマンフィルターの性能