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

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

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

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

load vdp_data.mat y

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

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

Реализуйте алгоритм сигма-точечного фильтра Калмана, чтобы оценить состояния осциллятора при помощи команд predict и correct. Вы сначала исправляете $$\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.

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

for k = 1:size(y)
    [CorrectedState,CorrectedStateCovariance] = correct(obj,y(k));
    [PredictedState,PredictedStateCovariance] = predict(obj);
end

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

В этом примере вы использовали 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, смотрите Используя, предсказывают и исправляют Команды.

Загрузите данные об ОДУ Ван дер Поля и задайте шаг расчета.

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

addpath(fullfile(matlabroot,'examples','ident','main')) % add example data

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')

rmpath(fullfile(matlabroot,'examples','ident','main')) % remove example data

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

Можно теперь оценить состояние системы с помощью команд correct и predict. Вы передаете значения 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)