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

Динамика объекта

Рассмотрите дискретный объект с аддитивным Гауссовым шумом $$w\left[n\right]$$ на входе $$u\left[n\right]$$:

Далее, позволить $$y_v\left[n\right]$$ будьте шумным измерением вывода $$y\left[n\right]$$, с $$v\left[n\right]$$ обозначение шума измерения:

Следующие матрицы представляют динамику объекта.

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];

Дискретный фильтр Калмана

Уравнения установившегося Фильтра Калмана для этой проблемы даны можно следующим образом.

В этих уравнениях:

Учитывая текущую оценку $$\underset{}{\overset{ˆ}{x}}\left[n|n\right]$$, обновление времени предсказывает значение состояния на следующей выборке n + 1 (предиктор "один шаг вперед"). Обновление измерения затем настраивает этот прогноз на основе нового измерения $$y_v[n+1]$$. Срок исправления является функцией инноваций, то есть, несоответствия между измеренными и ожидаемыми значениями $$y[n+1]$$. Этим несоответствием дают:

Инновационное усиление M выбрано, чтобы минимизировать установившуюся ковариацию ошибки оценки, учитывая шумовые ковариации:

Можно объединить время и уравнения обновления измерения в одну модель в пространстве состояний, Фильтр Калмана:

Этот фильтр генерирует оптимальную оценку $$\underset{}{\overset{ˆ}{y}}\left[n|n\right]$$ из $$y_n$$. Обратите внимание на то, что состояние фильтра $$\underset{}{\overset{ˆ}{x}}\left[n|n-1\right]$$.

Установившийся проект

Можно разработать установившийся Фильтр Калмана, описанный выше с функциональным kalman. Сначала задайте модель объекта управления с шумом процесса:

Здесь, первое выражение является уравнением состояния, и вторым является уравнение измерения.

Следующая команда задает эту модель объекта управления. Шаг расчета установлен в-1, чтобы отметить модель как дискретная, не задавая шаг расчета.

Plant = ss(A,[B B],C,0,-1,'inputname',{'u' 'w'},'outputname','y');

При предположении, что Q = R = 1, разработайте дискретный Фильтр Калмана.

Q = 1; 
R = 1;
[kalmf,L,P,M] = kalman(Plant,Q,R);

Эта команда возвращает модель в пространстве состояний kalmf фильтра, а также инновации получают M.

M

M = 3×1

    0.3798
    0.0817
   -0.2570

Входные параметры kalmf - вы и $$y_v$$, и. Выходные параметры являются объектом вывод и оценки состояния, $$y_e=\underset{}{\overset{ˆ}{y}}\left[n|n\right]$$ и $$\underset{}{\overset{ˆ}{x}}\left[n|n\right]$$.

Поскольку вы интересуетесь выходной оценкой $$y_e$$, выберите первый вывод kalmf и отбросьте остальных.

kalmf = kalmf(1,:);

Чтобы видеть, как фильтр работает, сгенерируйте некоторые входные данные и случайный шум и сравните отфильтрованный ответ $$y_e$$ с истинным ответом y. Можно или сгенерировать каждый ответ отдельно или сгенерировать обоих вместе. Чтобы моделировать каждый ответ отдельно, используйте lsim с одним только объектом первый, и затем с объектом и отфильтруйте присоединенный вместе. Объединенная альтернатива симуляции детализирована затем.

Блок-схема ниже показов, как сгенерировать и истинные и отфильтрованные выходные параметры.

Можно создать модель в пространстве состояний этой блок-схемы с функциями parallel и feedback. Сначала создайте модель комплектного оборудования с u, w, v как входные параметры и y и $$y_v$$ (измерения) как выходные параметры.

a = A;
b = [B B 0*B];
c = [C;C];
d = [0 0 0;0 0 1];
P = ss(a,b,c,d,-1,'inputname',{'u' 'w' 'v'},'outputname',{'y' 'yv'});

Затем используйте parallel, чтобы сформировать параллельную связь следующего рисунка.

sys = parallel(P,kalmf,1,1,[],[]);

Наконец, замкните круг датчика путем соединения объекта вывод $$y_v$$ отфильтровать вход $$y_v$$ с позитивными откликами.

SimModel = feedback(sys,1,4,2,1);   % Close loop around input #4 and output #2
SimModel = SimModel([1 3],[1 2 3]); % Delete yv from I/O list

Получившаяся имитационная модель имеет w, v, u как входные параметры и y и $$y_e$$ как выходные параметры. Просмотрите свойства InputName и OutputName проверить.

SimModel.InputName

ans = 3x1 cell array
    {'w'}
    {'v'}
    {'u'}

SimModel.OutputName

ans = 2x1 cell array
    {'y'  }
    {'y_e'}

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

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

n = length(t);
rng default
w = sqrt(Q)*randn(n,1);
v = sqrt(R)*randn(n,1);

Моделируйте ответы.

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

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

Сравните истинные и отфильтрованные ответы графически.

subplot(211), plot(t,y,'--',t,ye,'-'), 
xlabel('No. of samples'), ylabel('Output')
title('Kalman filter response')
subplot(212), plot(t,y-yv,'-.',t,y-ye,'-'),
xlabel('No. of samples'), ylabel('Error')

Первый график показывает истинный ответ y (пунктирная линия) и отфильтрованный вывод $$y_e$$ (сплошная линия). Второй график сравнивает погрешность измерения (штрихпунктир) с ошибкой оценки (тело). Этот график показывает, что уровень шума значительно уменьшался. Это подтверждено путем вычисления ошибок ковариации. Ошибочная ковариация прежде, чем отфильтровать (погрешность измерения):

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

MeasErrCov = 0.9992

Ошибочная ковариация после фильтрации (ошибка оценки) уменьшается:

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

EstErrCov = 0.4944

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

Изменяющийся во времени Фильтр Калмана является обобщением установившегося фильтра для изменяющихся во времени систем или систем LTI с неустановившейся шумовой ковариацией.

Рассмотрите следующее состояние объекта и уравнения измерения.

Изменяющийся во времени Фильтр Калмана дан следующими рекурсиями:

Здесь, $$\underset{}{\overset{ˆ}{x}}\left[n|n-1\right]$$ и $$\underset{}{\overset{ˆ}{x}}\left[n|n\right]$$ как описаны ранее. Дополнительно:

Для простоты были пропущены индексы, указывающие на временную зависимость матриц пространства состояний.

Учитывая начальные условия $$x\left[1|0\right]$$ и $$P\left[1|0\right]$$, можно выполнить итерации этих уравнений, чтобы выполнить фильтрацию. Необходимо обновить обоих оценки состояния $$x\left[n|.\right]$$ и ошибочные ковариационные матрицы $$P\left[n|.\right]$$ в каждый раз выборка.

Изменяющийся во времени проект

Чтобы реализовать эти рекурсии фильтра, сначала сгенерируйте шумные выходные измерения. Используйте шум процесса w и шум измерения v, сгенерированный ранее.

sys = ss(A,B,C,0,-1);
y = lsim(sys,u+w);      
yv = y + 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); 

for i = 1:length(t)
  % Measurement update
  Mn = P*C'/(C*P*C'+R);
  x = x + Mn*(yv(i)-C*x);   % x[n|n]
  P = (eye(3)-Mn*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,y,'--',t,ye,'-')
title('Time-varying Kalman filter response')
xlabel('No. of samples'), ylabel('Output')
subplot(212), plot(t,y-yv,'-.',t,y-ye,'-')
xlabel('No. of samples'), ylabel('Output')

Первый график показывает истинный ответ y (пунктирная линия) и отфильтрованный ответ $$y_e$$ (сплошная линия). Второй график сравнивает погрешность измерения (штрихпунктир) с ошибкой оценки (тело).

Изменяющийся во времени фильтр также оценивает ковариацию errcov ошибки оценки $$y-y_e$$ на каждой выборке. Постройте его, чтобы видеть, достиг ли ваш фильтр устойчивого состояния (как вы ожидаете со стационарным входным шумом).

subplot(211)
plot(t,errcov), ylabel('Error covar')

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

Сравните с ошибочной ковариацией оценки, выведенной от экспериментальных данных:

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

EstErrCov = 0.4934

Это значение меньше, чем теоретическое значение errcov и близко к значению, полученному для установившегося проекта.

Наконец, обратите внимание что окончательное значение $$M\left[n\right]$$ и установившееся значение M инновационной матрицы усиления совпадает.

Mn

Mn = 3×1

    0.3798
    0.0817
   -0.2570

M

M = 3×1

    0.3798
    0.0817
   -0.2570

Библиография

[1] Grimble, M.J., Устойчивое Управление производственным процессом: Оптимальный Подход Проекта для Полиномиальных Систем, Prentice Hall, 1994, p. 261 и стр 443-456.