Параметр и оценка состояния в Simulink Используя блок фильтра частиц

Этот пример использует:

Этот пример демонстрирует использование блока Particle Filter в Control System Toolbox™. Задача оценки параметра передаточной функции дискретного времени переформулирована и рекурсивно решена как проблема оценки состояния.

Введение

Control System Toolbox имеет три блока Simulink для нелинейной оценки состояния:

Фильтр частиц: Реализует алгоритм фильтра частиц дискретного времени.
Расширенный Фильтр Калмана: Реализует расширенный алгоритм Фильтра Калмана дискретного времени первого порядка.
Сигма-точечный фильтр Калмана: Реализует алгоритм сигма-точечного фильтра Калмана дискретного времени.

Эти блоки поддерживают оценку состояния с помощью нескольких датчиков, действующих на уровне различных частот дискретизации. Типичный рабочий процесс для использования этих блоков следующие:

Смоделируйте свой объект и поведение датчика с помощью функций Simulink или MATLAB.
Сконфигурируйте параметры блока.
Симулируйте фильтр и анализируйте результаты завоевать доверие в эффективности фильтра.
Разверните фильтр на своем оборудовании. Можно сгенерировать код для этих фильтров с помощью программного обеспечения Simulink Coder™.

Этот пример использует блок Particle Filter, чтобы продемонстрировать первые два шага этого рабочего процесса. Последние два шага кратко обсуждены в разделе Next Steps. Цель в этом примере состоит в том, чтобы оценить параметры передаточной функции дискретного времени (модель ошибки на выходе) рекурсивно, где параметры модели обновляются на каждом временном шаге, когда новая информация прибывает.

Если вы интересуетесь Расширенным Фильтром Калмана, смотрите пример "Оценочные состояния Нелинейной Системы с Несколькими, Многоскоростные Датчики". Использование Сигма-точечного фильтра Калмана выполняет подобные шаги к Расширенному Фильтру Калмана.

Добавьте папку в качестве примера в путь MATLAB.

addpath(fullfile(matlabroot,'examples','control','main'));

Моделирование объекта

Большинство алгоритмов оценки состояния использует модель объекта управления (функция изменения состояния), которые описывают эволюцию состояний объекта от одного временного шага до следующего. Эта функция обычно обозначается как $x[k+1]=f(x[k],w[k],u[k])$ , где x является состояниями, w является шумом процесса, u является дополнительными дополнительными входными параметрами, например, системные входные параметры или параметры. Блок фильтра частиц требует, чтобы вы обеспечили эту функцию в немного отличающемся синтаксисе $X[k+1]=f_{pf}(X[k],u[k])$ . Различия:

Фильтр частиц работает следующим траектории многих гипотез состояния (частицы), и блок передает все гипотезы состояния вашей функции целиком. Конкретно, если ваш вектор состояния имеет элементы, и вы выбираете частицы, чтобы использовать, имеет размерности $[N_s \; N_p]$ , где каждый столбец является гипотезой состояния.
Вы вычисляете удар шума процесса на гипотезах состояния в вашей функции $f_{pf}(...)$ . Блок не делает предположений на вероятностном распределении шума процесса и не нуждается как вход.

Функцией $f_{pf}(...)$ может быть функция MATLAB, которые выполняют ограничения MATLAB Coder™ или Функциональный блок Simulink. После того, как вы создадите $f_{pf}(...)$ , вы задаете имя функции в блоке Particle Filter.

В этом примере вы переформулируете проблему оценки параметра передаточной функции дискретного времени как проблему оценки состояния. Эта передаточная функция может представлять динамику процесса дискретного времени, или это может представлять некоторую динамику непрерывного времени вместе с сигналом reconstructor, таким как нулевой порядок, содержат. Примите, что вы интересуетесь оценкой параметров передаточной функции дискретного времени первого порядка:

$y[k] = \frac{20q^{-1}}{1-0.7q^{-1}} u[k] + e[k]$

Вот $y[k]$ является объект выходом, $u[k]$ вход объекта, $e[k]$ шум измерения, $q^{-1}$ оператор с временной задержкой так, чтобы $q^{-1}u[k]=u[t-1]$ . Параметризуйте передаточную функцию как $\frac{n\;q^{-1}}{1+d\;q^{-1}}$ , где $n$ и $d$ параметры, которые будут оценены. Передаточная функция и параметры могут быть представлены в необходимой форме пространства состояний несколькими способами, по выбору вектора состояния. Один выбор состоит в том $x[k]=[\; y[k]; d[k]; n[k] \;]$ , где вторые и третьи состояния представляют оценки параметра. Затем передаточная функция может быть эквивалентно записана как

$x[k+1] = \left[
 \begin{array}{c}
 -x_2[k] x_1[k] + x_3[k] u[k] \\
 x_2[k] \\
 x_3[k] \\
 \end{array}
 \right]$

Шумовая часть измерения $e[k]$ обработана в моделировании датчика. В этом примере вы реализуете выражение выше в функции MATLAB в векторизованной форме для вычислительного КПД:

type pfBlockStateTransitionFcnExample

function xNext = pfBlockStateTransitionFcnExample(x,u)
% pfBlockStateTransitionFcnExample State transition function for particle
%                                  filter, for estimating parameters of a 
%                                  SISO, first order, discrete-time transfer 
%                                  function model
%
% Inputs:
%    x - Particles, a NumberOfStates-by-NumberOfParticles matrix
%    u - System input, a scalar 
%
% Outputs:
%    xNext - Predicted particles, with the same dimensions as input x

% Implement the state-transition function
%    xNext = [x(1)*x(2) + x(3)*u; 
%             x(2); 
%             x(3)];
% in vectorized form (for all particles).
xNext = x;
xNext(1,:) = bsxfun(@times,x(1,:),-x(2,:)) + x(3,:)*u;

% Add a small process noise (relative to expected size of each state), to
% increase particle diversity
xNext = xNext + bsxfun(@times,[1; 1e-2; 1e-1],randn(size(xNext)));
end

Моделирование датчика

Блок Particle Filter требует, чтобы вы обеспечили функцию правдоподобия измерения, которая вычисляет вероятность (вероятность) каждой гипотезы состояния. Эта функция имеет форму $L[k] = h_{pf}(X[k],y[k],u[k])$ . $L[k]$ действительно ли $N_p$ элементом является вектор, где $N_p$ количество частиц, вы выбираете. $m^{th}$ элементом в $L[k]$ является вероятность $m^{th}$ частицы (столбец) в $X[k]$ . $y[k]$ измерение датчика. $u[k]$ дополнительный входной параметр, который может отличаться от входных параметров функции изменения состояния.

Датчик измеряет первое состояние в этом примере. Этот пример принимает, что ошибки между фактическими и предсказанными измерениями распределяются согласно Распределению Гаусса, но любое произвольное вероятностное распределение или некоторый другой метод могут использоваться, чтобы вычислить вероятности. Вы создаете $h_{pf}(...)$ , без обозначения даты задаете имя функции в блоке Particle Filter.

type pfBlockMeasurementLikelihoodFcnExample

function likelihood = pfBlockMeasurementLikelihoodFcnExample(particles, measurement)
% pfBlockMeasurementLikelihoodFcnExample Measurement likelihood function for particle filter
%
% The measurement is the first state
%
% Inputs:
%    particles   - A NumberOfStates-by-NumberOfParticles matrix
%    measurement - System output, a scalar
%
% Outputs:
%    likelihood - A vector with NumberOfParticles elements whose n-th
%                 element is the likelihood of the n-th particle

%#codegen

% Predicted measurement
yHat = particles(1,:);

% Calculate likelihood of each particle based on the error between actual
% and predicted measurement
%
% Assume error is distributed per multivariate normal distribution with
% mean value of zero, variance 1. Evaluate the corresponding probability
% density function
e = bsxfun(@minus, yHat, measurement(:)'); % Error
numberOfMeasurements = 1;
mu = 0; % Mean
Sigma = eye(numberOfMeasurements); % Variance
measurementErrorProd = dot((e-mu), Sigma \ (e-mu), 1);
c = 1/sqrt((2*pi)^numberOfMeasurements * det(Sigma));
likelihood = c * exp(-0.5 * measurementErrorProd);
end

Отфильтруйте конструкцию

Сконфигурируйте блок Particle Filter для оценки. Вы задаете имена функции правдоподобия изменения состояния и измерения, количество частиц и начальное распределение этих частиц.

Во вкладке System Model диалогового окна блока задайте следующие параметры:

Изменение состояния

Задайте функцию изменения состояния, pfBlockStateTransitionFcnExample, в Функции. При вводе имя функции и нажимаете Apply, блок обнаруживает, что ваша функция имеет дополнительный вход и создает входной порт StateTransitionFcnInputs. Вы соединяете свой системный вход с этим портом.

Инициализация

Задайте 10000 в Количестве частиц. Более высокое количество частиц обычно соответствует лучшей оценке по увеличенной вычислительной стоимости.
Задайте Gaussian в Распределении, чтобы получить начальный набор частиц от многомерного Распределения Гаусса. Затем задайте [0; 0; 0] в Среднем значении, потому что у вас есть три состояния и это - ваше лучшее предположение. Задайте diag([10 5 100]) в Ковариации, чтобы задать большое отклонение (неопределенность) в вашем предположении для третьего состояния и меньшем отклонении для первых двух. Очень важно, что этот начальный набор частиц распространен достаточно широкий (большое отклонение), чтобы покрыть потенциальное истинное состояние.

Измерение 1

Задайте имя своей функции правдоподобия измерения, pfBlockMeasurementLikelihoodFcnExample, в Функции.

Размер шага

В нижней части диалогового окна блока войдите 1 в Шаг расчета. Если у вас есть различный шаг расчета среди функций правдоподобия изменения состояния и измерения, или если у вас есть несколько датчиков с различными шагами расчета, они могут быть сконфигурированы в Block выходные параметры, вкладка Multirate.

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

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

Симуляция и результаты

Для простого теста истинная модель объекта управления симулирована с входными параметрами белого шума. Входные параметры и шумные измерения от объекта питаются блок Particle Filter. Следующая модель Simulink представляет эту настройку:

open_system('pfBlockExample')

Симулируйте систему и сравните предполагаемые и истинные параметры:

sim('pfBlockExample')
open_system('pfBlockExample/Parameter Estimates')

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

Поиск и устранение проблем

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

Поиск и устранение проблем фильтра частиц обычно выполняется путем изучения набора частиц и их весов, которые могут быть получены путем выбора Output все частицы и Выход все веса в Block выходные параметры, вкладка Multirate диалогового окна блока.

Первая проверка работоспособности должна гарантировать, что функции правдоподобия изменения состояния и измерения получают поведение вашей системы обоснованно хорошо. Если у вас есть имитационная модель для вашей системы (и следовательно доступ к истинному состоянию в симуляции), можно попытаться инициализировать фильтр истинным состоянием. Затем можно подтвердить, если функция изменения состояния вычисляет распространение времени истинного состояния точно и f, функция правдоподобия измерения вычисляет высокую вероятность для этих частиц.

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

Если точность оценки состояния прекрасна первоначально, но ухудшающийся в зависимости от времени, проблемой может быть степень вырождения частицы или обеднение частицы [2]. Степень вырождения частицы происходит, когда частицы распределяются слишком широко, в то время как обеднение частицы происходит из-за частиц, наносящих удар вместе после передискретизации. Степень вырождения частицы приводит к обеднению частицы в результате прямой передискретизации. Сложение искусственного шума процесса в функции изменения состояния, используемой в этом примере, является одним практическим подходом. Существует большое количество литературы по решению этих вопросов, и на основе вашего приложения более систематические подходы могут быть доступными. [1], [2] две ссылки, которые могут быть полезными.

Следующие шаги

Подтвердите оценку состояния: Если фильтр выполняет как ожидалось в симуляции, обычно эффективность далее подтверждена с помощью обширных симуляций Монте-Карло. Для получения дополнительной информации смотрите, Подтверждают Онлайновую Оценку состояния в Simulink. Можно использовать опции под группой Случайности в диалоговом окне блока Particle Filter, чтобы упростить эти симуляции.
Сгенерируйте код: блок Particle Filter поддерживает генерацию C и Кода С++ с помощью программного обеспечения Simulink Coder™. Функции, которые вы предоставляете этому блоку, должны выполнить ограничения программного обеспечения MATLAB Coder™ (если вы используете функции MATLAB, чтобы смоделировать вашу систему), и программное обеспечение Simulink Coder (если вы используете Функциональные блоки Simulink, чтобы смоделировать вашу систему).

Сводные данные

Этот пример показал, как использовать блок Particle Filter в Control System Toolbox. Вы оценили параметры передаточной функции дискретного времени рекурсивно, где параметры обновляются на каждом временном шаге, когда новая информация прибывает.

close_system('pfBlockExample', 0)

Удалите папку в качестве примера из пути MATLAB.

rmpath(fullfile(matlabroot,'examples','control','main'));

Ссылки

[1] Саймон, Дэн. Оптимальная оценка состояния: Кальман, H бесконечность и нелинейные подходы. John Wiley & Sons, 2006.

[2] Дусэ, Арно и Адам М. Йохансен. "Пример на фильтрации частицы и сглаживании: Пятнадцать лет спустя". Руководство нелинейной фильтрации 12.656-704 (2009): 3.

Документация