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

Открыть 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

Поддержка