В этом примере показано, как использовать расширенный Фильтр Калмана для обнаружения отказа. Пример использует расширенный Фильтр Калмана для онлайновой оценки трения простого двигателя постоянного тока. Существенные изменения в предполагаемом трении обнаруживаются и указывают на отказ. Этот пример использует функциональность от System Identification Toolbox™ и не требует Predictive Maintenance Toolbox™.
Двигатель моделируется как инерция J с коэффициентом демпфирования c, управляется крутящим моментом u. Моторная скорость вращения w и ускорение , измеренные выходные параметры.
Чтобы оценить коэффициент демпфирования c использование расширенного Фильтра Калмана, введите вспомогательное состояние для коэффициента демпфирования и обнулите его производную.
Таким образом, состояние модели, x = [w; c], и измерение, y, уравнения:
Уравнения непрерывного времени преобразовываются к дискретному времени с помощью приближения , где Ts является дискретным периодом выборки. Это дает уравнения модели дискретного времени, которые реализованы в pdmMotorModelStateFcn.m
и pdmMotorModelMeasurementFcn.m
функции.
Задайте параметры двигателя.
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. Шумовые условия являются дополнением.
Шум процесса и измерения имеет нулевое среднее значение, 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
:
Якобиан состояния задан в 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
Этот пример показал, как использовать расширенный Фильтр Калмана, чтобы оценить трение в простом двигателе постоянного тока и использовать оценку трения для обнаружения отказа.