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

Скрипт Open Live Script

В этом примере показано, как использовать расширенный Фильтр Калмана для обнаружения отказа. Пример использует расширенный Фильтр Калмана для онлайновой оценки трения простого двигателя постоянного тока. Существенные изменения в предполагаемом трении обнаруживаются и указывают на отказ. Этот пример использует функциональность от 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')

Figure contains 2 axes. Axes 1 with title Motor Outputs contains 2 objects of type line. These objects represent Measured Angular Velocity, Measured Angular Acceleration. Axes 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. Axes 1 with title Motor State - Velocity contains 3 objects of type line. These objects represent True value, Estimated value, Confidence interval. Axes 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. Axes 1 with title Velocity State Error contains an object of type line. Axes 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. Axes 1 with title Velocity Measurement Error contains an object of type line. Axes 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. Axes 1 with title Velocity Error Kurtosis contains an object of type line. Axes 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. The axes with title Friction Change Detection contains 2 objects of type line. These objects represent Estimated Friction, No-Fault Friction Bounds.

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

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

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

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

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

Создание расширенного фильтра Калмана

Симуляция

Расширенная эффективность фильтра Калмана

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

Похожие темы

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

Поддержка