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

Создайте расширенный объект фильтра Калмана для генератора. Используйте ранее написанные и сохраненные функции перехода и измерения, 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')

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

Загрузите данные Ван дер Поль и укажите шаг расчета.

vdpODEdata.mat содержит симуляцию Ван дер Поль с параметром нелинейности mu = 1, с использованием ode45, с начальными условиями [2;0]. Истинное состояние было извлечено со шаг расчета dt = 0.05.

load ('vdpODEdata.mat','xTrue','dt')
tSpan = 0:dt:5;

Получите измерения. В данном примере датчик измеряет первое состояние Гауссовым шумом со стандартным отклонением 0.04.

sqrtR = 0.04;
yMeas = xTrue(:,1) + sqrtR*randn(numel(tSpan),1);

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

myPF = particleFilter(@vdpParticleFilterStateFcn,@vdpMeasurementLikelihoodFcn);

Инициализируйте фильтр частиц в состоянии [2; 0] с единичной ковариацией и использованием 1000 частицы.

initialize(myPF,1000,[2;0],eye(2));

Выберите mean оценка состояния и systematic методы повторной дискретизации.

myPF.StateEstimationMethod = 'mean';
myPF.ResamplingMethod = 'systematic';

Оцените состояния, используя correct и predict команды и сохраните предполагаемые состояния.

xEst = zeros(size(xTrue));
for k=1:size(xTrue,1)
    xEst(k,:) = correct(myPF,yMeas(k));
    predict(myPF);
end

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

figure(1)
plot(xTrue(:,1),xTrue(:,2),'x',xEst(:,1),xEst(:,2),'ro')
legend('True','Estimated')

Рассмотрим нелинейную систему с входными 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);