Оцените состояния осциллятора Ван дер Поля с помощью расширенного алгоритма фильтра Калмана и измеренных выходных данных. Генератор имеет два состояния и один выход.

Создайте расширенный объект фильтра Калмана для генератора. Используйте ранее написанные и сохраненные функции перехода и измерения, vdpStateFcn.m и vdpMeasurementFcn.m. Эти функции описывают дискретное приближение к осциллятору Ван дер Поля с параметром нелинейности $\mathrm{mu}$ равным 1. Функции предполагают аддитивный процесс и измерительный шум в системе. Задайте начальные значения для двух состояний как [1; 0]. Это предположение для значения состояния в начальный момент времени k, основанный на знании выходов системы до времени k-1, $$\underset{}{\overset{ˆ}{x}}[k|k-1]$$.

obj = extendedKalmanFilter(@vdpStateFcn,@vdpMeasurementFcn,[1;0]);

Загрузите измеренные выходные данные y от генератора. В этом примере используйте моделируемые статические данные для рисунка. Данные хранятся в vdp_data.mat файл.

load vdp_data.mat y

Задайте технологические шумовые и измерительные шумовые ковариации генератора.

obj.ProcessNoise = 0.01;
obj.MeasurementNoise = 0.16;

Инициализируйте массивы, чтобы захватить результаты оценки.

residBuf = [];
xcorBuf = [];
xpredBuf = [];

Реализуйте расширенный алгоритм фильтра Калмана, чтобы оценить состояния генератора при помощи correct и predict команды. Вы сначала исправите $$\underset{}{\overset{ˆ}{x}}[k|k-1]$$ использование измерений во времени k добираться $$\underset{}{\overset{ˆ}{x}}[k|k]$$. Затем вы прогнозируете значение состояния на следующем временном шаге $$\underset{}{\overset{ˆ}{x}}[k+1|k]$$ использование $$\underset{}{\overset{ˆ}{x}}[k|k]$$оценка состояния на временном шаге k который оценивается с помощью измерений до k времени.

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

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

for k = 1:size(y)
    [Residual,ResidualCovariance] = residual(obj,y(k));
    [CorrectedState,CorrectedStateCovariance] = correct(obj,y(k));
    [PredictedState,PredictedStateCovariance] = predict(obj);
    
     residBuf(k,:) = Residual;
     xcorBuf(k,:) = CorrectedState';
     xpredBuf(k,:) = PredictedState';
   
end

Когда вы используете correct команда, obj.State и obj.StateCovariance обновляются с учетом скорректированного состояния и ошибки расчета состояния ковариации значений для временного шага k, CorrectedState и CorrectedStateCovariance. Когда вы используете predict команда, obj.State и obj.StateCovariance обновляются предсказанными значениями для временного шага k+1, PredictedState и PredictedStateCovariance. Когда вы используете residual команда, вы не изменяете никакие obj свойства.

В этом примере вы использовали correct перед predict потому что начальное значение состояния было $$\underset{}{\overset{ˆ}{x}}[k|k-1]$$, предположение для значения состояния в начальное время k на основе системных выходов до временного k-1. Если ваше начальное значение состояния $$\underset{}{\overset{ˆ}{x}}[k-1|k-1]$$, значение в предыдущий момент времени k-1 на основе измерений до k-1, затем используйте predict сначала команда. Для получения дополнительной информации о порядке использования predict и correct, см. Использование предсказания и исправление команд.

Постройте график предполагаемых состояний с помощью значений после посткоррекции.

plot(xcorBuf(:,1), xcorBuf(:,2))
title('Estimated States')

Постройте график фактического измерения, скорректированного оценочного измерения и невязки. Для функции измерения в vdpMeasurementFcn, измерение является первым состоянием.

M = [y,xcorBuf(:,1),residBuf];
plot(M)
grid on
title('Actual and Estimated Measurements, Residual')
legend('Measured','Estimated','Residual')

Оценка внимательно отслеживает измерение. После начального переходного процесса невязка остается относительно небольшим в течение всего запуска.

Рассмотрим нелинейную систему с входными u чье состояние x и измерительные y развиваться согласно следующим переходным и измерительным уравнениям:

Шум процесса w системы является аддитивным, в то время как шум измерения v является неаддитивным.

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

f = @(x,u)(sqrt(x+u));
h = @(x,v,u)(x+2*u+v^2);

f и h являются указателями на функцию анонимных функций, которые хранят функции перехода и измерения состояния, соответственно. В функции измерения, потому что шум измерения неаддитивен, v задается также как вход. Обратите внимание, что v задается как вход перед дополнительным входом u.

Создайте расширенный объект фильтра Калмана для оценки состояния нелинейной системы с помощью заданных функций. Задайте начальное значение состояния как 1, и шум измерения как неаддитивный.

obj = extendedKalmanFilter(f,h,1,'HasAdditiveMeasurementNoise',false);

Задайте ковариацию шума измерения.

obj.MeasurementNoise = 0.01;

Теперь можно оценить состояние системы с помощью predict и correct команды. Вы передаете значения u на predict и correct, которые в свою очередь передают их в функции переходов состояний и измерений, соответственно.

Исправьте оценку состояния с помощью y измерения[k] = 0.8 и входные u[k] = 0.2 в временной шаг k.

correct(obj,0.8,0.2)

Спрогнозируйте состояние на следующем временном шаге, задайте u[k] = 0.2.

predict(obj,0.2)

Найдите ошибку, или невязку, между предсказанием и измерением.

[Residual, ResidualCovariance] = residual(obj,0.8,0.2);