Обнаружение отказа Используя расширенный фильтр Калмана

В этом примере показано, как использовать расширенный Фильтр Калмана в обнаружении отказа. Пример использует расширенный Фильтр Калмана в онлайновой оценке трения простого двигателя постоянного тока. Существенные изменения в предполагаемом трении обнаруживаются и указывают на отказ. Этот пример использует функциональность от System Identification Toolbox™ и не требует Predictive Maintenance Toolbox™.

Модель электродвигателя

Двигатель моделируется как инерция J с затуханием коэффициента c, управляется крутящим моментом u. Моторная скорость вращения 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

Симулируйте модель с помощью уравнения обновления состояния и добавьте шум процесса в состояния модели. Десять секунд в симуляцию, обеспечьте изменение в моторном трении. Используйте функцию измерения модели, чтобы симулировать моторные датчики и добавить шум измерения в выходные параметры модели.

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

Ответы ввода - вывода модели указывают, что это затрудняет, чтобы обнаружить изменение трения непосредственно от измеренных сигналов. Расширенный Фильтр Калмана позволяет нам оценить состояния, в частности состояние трения. Сравните истинные состояния модели и оцененные состояния. Предполагаемые состояния показывают с доверительными интервалами, соответствующими 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, 
subplot(211),plot(t,xSigTrue(1,:)-xSigEst(1,:))
title('Velocity State Error')
subplot(212),plot(t,xSigTrue(2,:)-xSigEst(2,:))
title('Friction State Error')

Диаграммы погрешностей показывают, что фильтр адаптируется после изменения трения в 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')

Ускоряющая диаграмма погрешностей показывает незначительные различия в средней погрешности приблизительно 10 секунд, когда отказ введен. Просмотрите ошибочную статистику, чтобы видеть, может ли отказ быть обнаружен от вычисленных ошибок. Ускорение и ошибки скорости, как ожидают, будут нормально распределены (шумовые модели являются все Гауссовыми). Поэтому эксцесс ускоряющей ошибки может помочь идентифицировать, когда изменение распределения ошибок от симметричного до асимметричного из-за трения изменяется и получившееся изменение в распределении ошибок.

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

При игнорировании первых 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

Сводные данные

Этот пример показал, как использовать расширенный Фильтр Калмана, чтобы оценить трение в простом двигателе постоянного тока и использовать оценку трения в обнаружении отказа.

Документация