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

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

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

Загрузите данные ODE van der Pol и укажите время выборки.

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

addpath(fullfile(matlabroot,'examples','control','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','control','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, которые, в свою очередь, передают их в функции перехода и измерения состояния соответственно.

Корректировка оценки состояния с помощью измерения 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);