Обнаружение неисправностей с помощью расширенного фильтра Калмана

Открыть сценарий в реальном времени

В этом примере показано, как использовать расширенный фильтр Калмана для обнаружения неисправностей. В примере используется расширенный фильтр Калмана для оперативной оценки трения простого двигателя постоянного тока. Обнаруживаются значительные изменения в расчетном трении и указывают на неисправность. В этом примере используются функциональные возможности Toolbox™ идентификации системы и не требуется Toolbox™ предиктивного обслуживания.

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

Двигатель моделируется как инерция J с коэффициентом демпфирования c, приводимым в действие крутящим моментом u. Угловая скорость w двигателя и ускорение $_{}^{w˙}$ являются измеренными выходами.

$_{}^{} w˙= (u-cw)$ /J

Чтобы оценить коэффициент демпфирования c с помощью расширенного фильтра Калмана, введите вспомогательное состояние для коэффициента демпфирования и установите его производную в нуль.

$_{}^{} c˙=0$

Таким образом, состояние модели, x = [w; c], и измерения, y, уравнения:

$[\begin{array}{c} _{}^{} \\ _{}^{} \end{array} w˙c˙]= \begin{array}{c} [(u-cw \\ )/ \end{array}$ J0]

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

Уравнения непрерывного времени преобразуются в дискретное время с использованием аппроксимационного $_{}^{} \frac{_{}_{}}{_{x˙=xn+1-xnTs}}$ , где Ts - период дискретной выборки. Это дает уравнения модели дискретного времени, которые реализованы в pdmMotorModelStateFcn.m и pdmMotorModelMeasurementFcn.m функции.

$[\begin{array}{c} _{wn +} \\ _{1cn} \end{array} + \begin{array}{c} 1_{]} =_{[} {wn}_{} +_{} (_{} \\ _{} \end{array} un-cnwn$ ) Ts/Jcn]

$_{yn} = \begin{array}{c} _{[} \\ {wn}_{} (_{}_{} un-cnwn \end{array}$ )/J]

Укажите параметры двигателя.

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} _{wn +} \\ _{1cn} \end{array} + \begin{array}{c} 1_{]} =_{[} {wn}_{} +_{} (_{} \\ _{} \end{array} un-cnwn)$ Ts/Jcn] + q

$_{yn} = \begin{array}{c} _{[} \\ {wn}_{} (_{}_{} un-cnwn \end{array}) /$ J] + 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{}{} \begin{array}{cc} _{}_{} & _{}_{} \\ ∂∂xf=[1-Tscn/J-Tswn/J01 \end{array}]$

$\frac{}{} \begin{array}{cc} _{} & _{} ∂∂xh=[10-cn/J-wn/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.

Резюме

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

Связанные темы

Индикаторы условий на основе модели

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