Стационарный фильтр Калмана

Рассмотрим следующую дискретную установку с гауссовым шумом w на входе и измерительным шумом v на выходе:

Целью является проектирование фильтра Калмана для оценки истинного выхода установки ${\mathit{}}_{\mathit{yt}}\mathit{[}n]\mathit{}=\mathit{}y\mathit{[}n\mathit{]}$-v [n] на основе шумных $\mathit{}\mathit{}$измерений y [n]. Этот стационарный фильтр Калмана использует следующие уравнения для этой оценки.

Обновление времени:

Обновление измерений:

Здесь,

(Эти уравнения обновления описывают current устройство оценки типа. Для получения информации о различии между current оценщики и delayed оценщики, см. kalman.)

Проектирование фильтра

Вы можете использовать kalman функция для проектирования этого стационарного фильтра Калмана. Эта функция определяет оптимальный стационарный коэффициент усиления М фильтра для конкретной установки на основе ковариации Q технологического шума и ковариации R шума датчика, которую вы обеспечиваете. В этом примере используются следующие значения для матриц состояния и пространства завода.

A = [1.1269   -0.4940    0.1129 
     1.0000         0         0 
          0    1.0000         0];

B = [-0.3832
      0.5919
      0.5191];

C = [1 0 0];

D = 0;

Для этого примера задайте G = B, что означает, что шум w процесса является входным шумом присадки. Кроме того, установка H = 0 означает, что входной шум w не оказывает прямого влияния на выходной сигнал y. Эти допущения дают более простую модель установки:

Когда H = 0, можно показать, что ${\mathit{}}_{\mathit{My}}\mathit{=}{\mathit{}}_{\mathit{}}$CMx ( см.kalman). Вместе эти предположения также упрощают уравнения обновления для фильтра Калмана.

Обновление времени:

Обновление измерений:

Для проектирования этого фильтра сначала создайте модель завода с вводом для w. Установите время выборки равным -1 пометить установку как дискретную (без определенного времени выборки).

Ts = -1;
sys = ss(A,[B B],C,D,Ts,'InputName',{'u' 'w'},'OutputName','y');  % Plant dynamics and additive input noise w

Ковариация технологического шума Q и ковариация шума датчика R значения больше нуля, которые обычно получаются из исследований или измерений системы. В этом примере задайте следующие значения.

Q = 2.3; 
R = 1; 

Используйте kalman для проектирования фильтра.

[kalmf,L,~,Mx,Z] = kalman(sys,Q,R);

Эта команда создает фильтр Калмана, kalmfмодель «состояние-пространство», которая реализует уравнения «время-обновление» и «измерение-обновление». Входами фильтра являются вход u установки и выход y шумной установки. Первый выход kalmf является оценочной $\stackrel{yˆ}{\mathit{}}$ истинного выхода установки, а остальные выходные данные являются оценками состояния $\stackrel{xˆ}{\mathit{}}$.

В этом примере отбрасывайте оценки состояния и сохраняйте только первый выход, $\stackrel{yˆ}{\mathit{}}$.

kalmf = kalmf(1,:);

Использовать фильтр

Чтобы увидеть, как работает этот фильтр, создайте некоторые данные и сравните отфильтрованный отклик с истинным откликом завода. Вся система показана на следующей схеме.

Для моделирования этой системы используйте sumblk для создания входных данных для измерения шума v. Затем используйте connect присоединиться sys и фильтр Калмана вместе, чтобы u является общим входом и шумным выходом установки y подается на другой вход фильтра. Результатом является расчетная модель с входными данными. w, v, и u и результаты yt (истинный ответ) и ye (отфильтрованный или оцененный $\stackrel{yˆ\; ответа}{\mathit{}}$). Сигналы yt и ye являются выходами установки и фильтра соответственно.

sys.InputName = {'u','w'};
sys.OutputName = {'yt'};
vIn = sumblk('y=yt+v');

kalmf.InputName = {'u','y'};
kalmf.OutputName = 'ye';

SimModel = connect(sys,vIn,kalmf,{'u','w','v'},{'yt','ye'});

Чтобы смоделировать поведение фильтра, создайте известный синусоидальный входной вектор.

t = (0:100)';
u = sin(t/5);

Создание векторов шума технологического процесса и шума датчика с использованием одинаковых значений ковариации шума Q и R который использовался при проектировании фильтра.

rng(10,'twister');
w = sqrt(Q)*randn(length(t),1);
v = sqrt(R)*randn(length(t),1);

Наконец, смоделировать ответ с помощью lsim.

out = lsim(SimModel,[u,w,v]);

lsim генерирует отклик на выходах yt и к входам, применяемым при w, v, и u. Извлеките yt и вы канализируете и вычисляете измеренный отклик.

yt = out(:,1);   % true response
ye = out(:,2);  % filtered response
y = yt + v;     % measured response

Сравните истинный ответ с отфильтрованным.

clf
subplot(211), plot(t,yt,'b',t,ye,'r--'), 
xlabel('Number of Samples'), ylabel('Output')
title('Kalman Filter Response')
legend('True','Filtered')
subplot(212), plot(t,yt-y,'g',t,yt-ye,'r--'),
xlabel('Number of Samples'), ylabel('Error')
legend('True - measured','True - filtered')

Как показывает второй график, фильтр Калмана уменьшает ошибку yt - y за счет измерения шума. Чтобы подтвердить это уменьшение, вычислите ковариацию ошибки перед фильтрацией (ковариация ошибки измерения) и после фильтрации (ковариация ошибки оценки).

MeasErr = yt-yt;
MeasErrCov = sum(MeasErr.*MeasErr)/length(MeasErr)

MeasErrCov = 0

EstErr = yt-ye;
EstErrCov = sum(EstErr.*EstErr)/length(EstErr)

EstErrCov = 0.3479

Изменяющаяся во времени конструкция фильтра Калмана

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

Изменяющийся во времени фильтр Калмана имеет следующие уравнения обновления. В изменяющемся во времени фильтре и ковариация $\mathit{ошибок\; P}\mathit{[}n$], и инновационный коэффициент усиления ${\mathit{}}_{\mathit{}}Mx\mathit{}[$n] могут изменяться со временем. Уравнения обновления времени и измерений можно изменить для учета изменения времени следующим образом. ( см.kalman для получения более подробной информации об этих выражениях.)

Обновление времени:

Обновление измерений:

Можно реализовать изменяющийся во времени фильтр Калмана в Simulink ® с помощью блока фильтра Калмана. Пример, демонстрирующий использование этого блока, см. в разделе Оценка состояния с использованием изменяющегося во времени фильтра Калмана. Для этого примера примените изменяющийся во времени фильтр в MATLAB ®.

Чтобы создать изменяющийся во времени фильтр Калмана, сначала создайте шумную реакцию установки. Моделирование реакции установки на входной сигнал u и шум технологического процесса w определено ранее. Затем добавьте шум измерения v на смоделированный истинный ответ yt для получения шумового ответа y. В этом примере ковариации векторов шума w и v не меняйтесь со временем. Однако можно использовать ту же процедуру для нестационарного шума.

yt = lsim(sys,[u w]);   
y = yt + v;         

Затем реализуйте уравнения обновления рекурсивного фильтра в for цикл.

P = B*Q*B';         % Initial error covariance
x = zeros(3,1);     % Initial condition on the state

ye = zeros(length(t),1);
ycov = zeros(length(t),1); 
errcov = zeros(length(t),1); 

for i=1:length(t)
  % Measurement update
  Mxn = P*C'/(C*P*C'+R);
  x = x + Mxn*(y(i)-C*x);   % x[n|n]
  P = (eye(3)-Mxn*C)*P;     % P[n|n]

  ye(i) = C*x;
  errcov(i) = C*P*C';

  % Time update
  x = A*x + B*u(i);        % x[n+1|n]
  P = A*P*A' + B*Q*B';     % P[n+1|n] 
end

Сравните истинный ответ с отфильтрованным.

subplot(211), plot(t,yt,'b',t,ye,'r--')
xlabel('Number of Samples'), ylabel('Output')
title('Response with Time-Varying Kalman Filter')
legend('True','Filtered')
subplot(212), plot(t,yt-y,'g',t,yt-ye,'r--'),
xlabel('Number of Samples'), ylabel('Error')
legend('True - measured','True - filtered')

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

figure
plot(t,errcov) 
xlabel('Number of Samples'), ylabel('Error Covariance'),

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

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

MeasErr = yt - y;
MeasErrCov = sum(MeasErr.*MeasErr)/length(MeasErr)

MeasErrCov = 0.9871

EstErr = yt - ye;
EstErrCov = sum(EstErr.*EstErr)/length(EstErr)

EstErrCov = 0.3479

Наконец, когда изменяющийся во времени фильтр достигает устойчивого состояния, значения в матрице усиления Mxn соответствуют вычисленным kalman для установившегося фильтра.

Mx,Mxn

Mx = 3×1

    0.5345
    0.0101
   -0.4776

Mxn = 3×1

    0.5345
    0.0101
   -0.4776