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

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

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;

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

Откорректируйте оценку состояния измерением yK]=0.8 и вход uK]=0.2 на временном шаге k.

correct(obj,0.8,0.2)

Предскажите состояние на следующем временном шаге, учитывая uK]=0.2.

predict(obj,0.2)

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

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