Оптимальное оценивание случайного процесса по критерию максимума апостериорной вероятности Текст научной статьи по специальности «Математика»

УДК 629.7+531.383

А. А. ЛОБАТЫЙ, Ю. Ф. ЯЦЫНА, Н. Н. АРЕФЬЕВ

ОПТИМАЛЬНОЕ ОЦЕНИВАНИЕ СЛУЧАЙНОГО ПРОЦЕССА ПО КРИТЕРИЮ МАКСИМУМА АПОСТЕРИОРНОЙ ВЕРОЯТНОСТИ

Белорусский национальный технический университет

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

Априорная плотность вероятности оцениваемого процесса изначально считается гауссовой дифференцируемой функцией, что позволяет разложить её в ряд Тейлора без использования в промежуточных преобразованиях характеристических функций и разложения на гармоники. Для малых интервалов времени плотность вероятности вектора ошибок измерений по определению так же задается гауссовой с нулевым математическим ожиданием. Это даёт возможность получить математическое выражение для функции невязки, характеризующей отклонение значений реального измерения процесса от его математической модели.

Для определения оптимальной апостериорной оценки вектора состояния задается предположение, что эта оценка соответствует ее математическому ожиданию - максимуму апостериорной плотности вероятности. Это даёт возможность на основе формулы Байеса для априорной и апостериорной плотности вероятности получить уравнение Стратоновича-Кушнера.

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

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

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

Эволюция вектора состояния Х(X) таких систем в общем случае может описываться нелинейным векторно-матричным уравнением вида

X (X) = ф( X ,и, X), X (¿0) = Хо, (1)

где X = [х1, ..., хп]Т = Х(х) - п-мерный вектор фазовых координат системы, ф - нелинейная векторная функция, Хо - случайный вектор начального состояния, и(X) - детерминированный вектор управлений (внешних воздействий).

Источники информации о состоянии системы представляют собой вектор 2 (X) размерности т < п . Математическая модель совокупности измерителей в общем случае может описываться стохастическим векторно-ма-тричным уравнением в форме Ланжевена

£(X) = №,X,X) + Н(X) .£(0, 2(Хо) = 2й. (2)

Здесь т - нелинейная векторная функция Н(X) - матрица коэффициентов, ^(0 - вектор случайных ошибок измерений. Большинство современных высокоточных измерителей можно считать безынерционными, а их инструмен-

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

Задача состоит в том, чтобы на основе вектора измерений Z(t) и известных математических моделей (1) и (2) получить оптимальную по заданному критерию оценку вектора состояния X (t) объекта управления или части фазовых координат вектора X(t), чтобы впоследствии использовать X (t ) для управления объектом.

В зависимости от вида функций и критерия оптимизации могут быть построены конкретные алгоритмы получения оценки X (t). Наибольшее распространение получили алгоритмы, основанные на линейной форме математической модели системы (1) и безынерционном измерителе (2). Такие математические модели объекта управления и измерителя имеют вид [3, 4]

X(t) = D(t)X(t) + U(t) + B(t)Ç(t), X(to ) = Xo, (3) Z (t ) = С (t ) • X (t ) + H (t ) -Z(t ). (4)

Здесь D(t), B(t), С (t), H (t) - матрицы коэффициентов, Ç(t), Z(t) - векторы некоррелированных белых гауссовых шумов с нулевыми математическими ожиданиями и матрицами интенсивностей G(t ) и N (t ), соответственно.

В качестве критерия оптимизации чаще всего рассматривается минимум среднего квадрата ошибки оценивания min M [( X (t ) - X (t ))2], где М[...] - символ математического ожидания.

Заметим, что практически любые математические модели систем с мультипликативными и небелыми шумы путем расширения вектора состояния и использования формирующих фильтров [5] можно свести к моделям с белыми аддитивными шумами. Модели вида (3) и (4) позволяют получить достаточно строгий вывод алгоритма вычисления оптимальной оценки X (t).

Модель измерителей (4) во многих случаях можно считать адекватной реальным устройствам (датчикам). В то же время линейная модель объекта управления (3) как правило, существенно отличается от нелинейной модели (1). Ошибки линеаризации, методические допущения порождают неопределенности в построении моделей процесса X(t).

Случайные процессы вида (3) с аддитивными белыми шумами являются марковскими

(точнее - марковскими первого порядка) [6], так как закон распределения такого процесса зависит только от его значений в текущий момент времени (непосредственно предшествующий будущему) и не зависит от того, какие значения принимал этот процесс в момент времени, предшествующий данному. Для марковского процесса плотность распределения вероятности (ПРВ) вектора фазовых координат Х(?) полностью определяется двумя функциями -первой функцией плотности вероятности /1(Х, ?) и плотностью вероятности перехода /1(Х(^1)|Х(^)). В некоторых случаях марковский процесс удобно описывать с помощью характеристических функций, связанных с плотностями У1(Х, ?) и /1(Х(^1) |Х(?)) через преобразование Фурье [7].

Непрерывный марковский процесс (3) может быть построен как предел суммы большого числа малых приращений. В соответствии с центральной предельной теоремой теории вероятностей закон распределения суммы индивидуально малых случайных величин с ростом числа слагаемых сводится к нормальному закону распределения [8]. Следовательно, для многомерного случайного марковского процесса Х(?) ПРВ имеет вид

/ (X, ?) = - 1

|0|

(5)

хехр{-1 (X(X) - тх (х))Т ©-1 (X(X) - тх (х))|.

В выражении (5) | © | - определитель ковариационной матрицы © вектора Х(х), тх (X) -математическое ожидание вектора Х(х). В скалярной форме выражение (5) имеет вид

/ (Х, 0 - 1

I© I

(6)

xexp< -

21 © I

Ё ©pq (xp - mxp )(Xq - mxq )

pq=1

xq '

где ©^ - алгебраическое дополнение элемента ©рд в определителе | © |.

Рассмотрим эволюцию ПРВ вектора Х(х) на бесконечно малом (элементарном) интервале времени Ах . Разложим функцию /(X, X) вида (6) в ряд Тейлора в окрестностях точки ? + Ах .

/ (х, ?) = / (Х, ? + АХ) +

п д

+ ^ [Атк (X, X) / (X, X + Ах] +

к=1 дхк

-[ЛипЫг (X, X) f (X, X нЛX ]н...

1 п д2

н- £ т-НЛти(X,X)f(X,Xнлх]н

2! к,/=1 дхк дх1 (7)

1 » д2

н— £ -

3! к ,1 ,г=1 дхк дх1 дхг

В выражении (7) через Ьт обозначены условные математические ожидания и условные корреляционные моменты приращений компонент хк (X) вектора X(X) при изменении аргумента X на малом интервале Ах.

Атк (X, X) = М [ хк (X + AX) - хк (X)] =

Ак (X, X )Ах + о(Ах),

Аты (X, X) =М [(хк (X + AX) - хк (X))(хг (X + AX) -

х1 (X ))] = Вы (X, X )А + o(Ax), (8 )

Атк1г (X, X) =М [(хк (X + AX) - хк (X))(х1 (X + AX) -

- х1 (X))(хг (X + Ах) - хг (X)] = о( Ах),

Атк1г...( X, X) = о(Ы),

где о(Ах) - величина высшего порядка малости чем Ах; к, I, г,... = 1,2,3,...

В выражениях (8) Ак (X, X) - компоненты вектора сноса А( X, X), Вк1 (X, X) - компоненты матрицы диффузии В( X, X) процесса X (X) вычисляются по формулам [9]

" хк (X + Ах) -х к (х)

Ак (X, X) = Нт М

Ах ^0

Ах

= Нт М

Ах ^0

Вк1 (X, X) = (хк (X+Ах )-хк (X))(хг (X+Ах )-х1 (X)) Ах

(9)

(10)

д/ (X, X) =_¥Т

дх

2 х

А( X, X) / (X, X)

УТ (ВТ (X, X) f (X, X))

(11)

где УТ =

5х1 дхп

- векторный оператор диф-

ференцирования.

Иногда уравнение (11) записывают в форме

/ (X, X)

дх

= - div п( X, X),

(12)

где = Ух , п(X, X) - плотность потока вероятности вида

п( X, X) = А( X, X) / (X, X) -

-1УТ 2 х

УТ (ВТ (X, X) f (X, X))

(13)

Уравнение Фоккера-Планка-Колмогорова (11) или (12) является уравнением в частных производных параболического типа. Его решение представляет значительную трудность, однако на основе его использования можно решить ряд задач, в частности - задачу апостериорного оценивания фазовых координат системы, описываемой марковской математической моделью.

Для определения апостериорной плотно-

л

сти вероятности /(X, X + Ах) процесса X(X), описываемого выражением (3) с учетом математической модели измерений (4) используем формулу Байеса [9] для плотностей вероятности

Условные вероятностные моменты Ат порядков выше второго имеют порядок малости о(Ах) в соответствии с определением [6] марковского процесса.

Разделив обе части выражения (7) на Ах и перейдя к пределу Ах ^ 0 с учетом выражений (9)-(10) и того, что

Ит / (X + Ах) -/ (X, X) = д/ (X, X) лх^0 Ах дх

получим многомерное уравнение Фоккера-План-ка-Колмогорова для априорной ПРВ процесса дх) [3, 5, 9].

/ (X, X + Ах) = / (X, X + Ах) / (2, X + Ах | X (X))

ТО

| /(X, X + Ах)/(2, X + Ах | X (X))dX

(14)

В выражении (14) /(X, X + Ах) - априорная ПРВ, эволюция которой в пространстве и во времени описывается уравнением (12)./(2, X + Ах | X(X)) - условная ПРВ вектора 2(X + Ах) при фиксированном состоянии вектора Х(х) в момент времени X.

Для малых интервалов времени Ах ПРВ входящего в выражение (4) вектора ошибок измерений ^(х + Ах) по определению [4] является гауссовой с нулевым математическим ожиданием. В соответствии с выражением (4) ПРВ /(2, X | X(X)) также является гауссовой с мате-

Л л

матическим ожиданием С (X) X (X), где X (X) -математическое ожидание апостериорной ПРВ

л л

/(X, X). Для симметричной ПРВ /(X, X) век-

—оо

тор X(X) соответствует максимуму /(X, X).

л

Следовательно, ПРВ /(2, X + А/ | X(X)) имеет вид

/(2 (X + АХ) | X (X)) =

1

^/(2гсГ |б(/ + АX)|

(15)

1

хехр ^ — Р( X, 2, X + АX Н,

2(X + АХ) - С(X + Q) X(X + АХ)

(16)

xQ _1(Х + АХ)

2 (X + АХ) - С (X + Q) X (X + АХ)

/ (X, X) = / (X, X) х

ехр - 2 Р( X, 2, X)

то 1 /(X, X )ехр —то — 2 р( X, 2, X) dX

(17)

1 + - 2 р( X, 2, X)

ТО | /(X, X) —ТО 1 — 2 р( X, 2, X) dX

/ (X, X + А/) = / (X, X) - А/ • сИУ п( X, X) -

1 Л Л <» Л

—АХр(X, г, X) f (X, X) + АХ Г div п(X, X)dX + (19)

2

—да

1 Л да Л Л

+-АXf(X, X) Г р(X, г, X) f (X, Х^,

2 —да

где Шу п( X, X) = п( X, X).

Предполагая существование предела

где |0(/ + A/)| - определитель дисперсионной матрицы процесса + ^) (интенсивность шума), р( X, 2, X + AX) - функция невязки [3, 4, 5], в данном случае вычисляемая по формуле

р( X, 2, X + АX) =

Нт —

AX^0 Ax

/ (X, X + AX) -/ (X, X)

д/ (X, X) дx

, (20)

после предельного перехода получим интегро-дифференциальное уравнение Стратоновича-Кушнера для апостериорной плотности вероятности [3, 5, 9]

С учетом (15) выражение (14) для момента X представим в виде

_ 1 2 где

д/(X,X) .

у 4 ' = _СТУ п(X, X) _ дX

Л да л л

р(X, 2, X) _ |р(X, 2, X) /(X, X)сК

п( X, X) = А( X, X) / (X, X) -

VI (ВТ (X, X) f (X, X)

(21)

/ (X, X),

(22)

Рассматривая эволюцию ПРВ /(X, X) на элементарном малом интервале Ax разложим экспоненциальную функцию по времени в ряд Маклорена (ех «1 + X + о(ф .

/ (X, X + A/) = / (X, X + A/) -

(18)

Использование уравнение (21) позволяет решать различные задачи [3, 4, 5, 9, 11] фильтрации, идентификации, прогноза состояния стохастических систем.

Умножив обе части выражения (21) на X (X) и проинтегрировав по X в бесконечных пределах, получим апостериорное уравнение для

л

математического ожидания М [ X (X)] = X (X), которое является оптимальной оценкой вектора X( X) по критерию максимума апостериорной вероятности.

X(X) = | А(X, X) /(X, X) -

В выражении (18) учтем, что на интервале A/ априорная ПРВ /(X, X+А/) =f(X,X)+A/ div п^,/), п^, X) - априорная плотность потока вероятности процесса X(x) вида (13).

После соответствующих преобразований выражения (18), удерживая члены первого порядка малости относительно А/, получим выражение [9]

1 да

-1 I

2

1

УХ (ВТ (X, X) f (X, X)

dX -

(23)

-- | (X(X) -X(X))р(X,2,X)/(X,/^.

2 -да

В выражении (23) учтено, что при X( X) = ±<х>,

/(X, X)=0 и етм=0.

д*

X

—то

то

В простейшем случае линейной постановки задачи, когда процесс X (X) описывается уравнением (3), а измеритель - выражением (4), вектор сноса имеет вид А( X, X) = В(х) X (X) + и (X), а матрица диффузии равна интенсивности шума £(х) В(X, X) = G(X). В этом случае выражение (23) принимает вид [3, 4]

X(X) = D(X) X(X) + и(X) + R(X)CТ (X^- (X) х (24)

х^(X) - С(X)X(X)],X= М[X0 ].

В выражении (24) Я (X) - апостериорная корреляционная матрица, дифференциальное уравнение для которой получается умножением

выражения (21) на [X(X) -X(х)]^(X) -X(х)]Т и интегрированием по X в бесконечных пределах.

Я(х) = В(х) Я(х) + Я(х) ВТ (X) + О (X) -

- R(X )СТ (X ^-1(x )С (X) R(X), R(Xо) = Rо.

(25)

Выражения (24) и (25) представляют собой классический фильтр Калмана (Калмана-Бью-си) [3, 4].

Заметим, что точно такой же результат получается при выводе фильтра Калмана (24)-(25) для линейной постановки задачи по критерию минимума среднего квадрата ошибки на основе использования интегрального уравнения Винера-Хопфа [5]. Однако использование дифференциального уравнения Стратоно-вича-Кушнера (21) даёт больше возможностей для получения различных алгоритмов апостериорной обработки информации [10, 11, 12].

ВЫВОД

Таким образом, получен оригинальный вывод уравнения Стратоновича-Кушнера для апо-

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

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

Для определения оптимальной апостериорной оценки вектора состояния задаем, что эта оценка соответствует ее математическому ожиданию - максимуму симметричной апостериорной плотности вероятности. Это даёт возможность на основе формулы Байеса для априорной и апостериорной плотности вероятности получить уравнение Стратоновича-Кушнера.

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

Литература

1. Красильщиков, М. Н. Современные информационные технологии в задачах навигации и наведения беспилотных маневренных летательных аппаратов / М. Н. Красильщиков, Г. Г. Серебряков. - М.: Физматлит, 2009. - 558 с.

2. Алешин, Б. С. Ориентация и навигация подвижных объектов: современные информационные технологии / Б. С. Алешин, К. К. Веремеенко, А. И. Черноморский. - М.: Физматлит, 2006. - 424 с.

3. Балакнишнан, А. В. Теория фильтрации Калмана: Пер. с англ. / А. В. Балакнишнан. - М.: Мир, 1988. - 1000 с.

4. Синицин, И. Н. Фильтры Калмана и Пугачева / И. Н. Синицин. - М.: Университетская книга, 2006. - 640 с.

5. Пупков, К. П. Методы классической и современной теории автоматического управления: в 5 т. / К. П. Пупков, Н. Д. Егупов. - М.: Издательство МГТУ им. Н. Э. Баумана, 2004. - 656 с.

6. Тихонов, В. И. Марковские процессы / В. И. Тихонов, И. Н. Миронов. - М. А.: Советское радио, 1977. - 450 с.

7. Пугачев, В. С. Теория стохастических систем / В. С. Пугачев, И. Н. Синицын. - М.: Логос, 2004. - 1000 с.

8. Вентцель, Е. С. Теория вероятностей / Е. С. Вентцель. - М.: Высшая школа, 1999. - 576 с.

9. Казаков, И. Е. Методы оптимизации стохастических систем / И. Е. Казаков, Д. И. Гладков. - М.: Наука, 1987. -304 с.

10. Лобатый, А. А. Особенности применения фильтров Калмана-Бьюси в комплексах ориентации и навигации / А. А. Лобатый, А. С. Бенкафо. - М.: Доклады БГУИР, 2013. - С. 67-71.

11. Лобатый, А. А. Оценка навигационных параметров подвижного объекта в условиях многорежимности / А. А. Лобатый, А. С. Бенкафо. - М.: Доклады БГУИР, 2014. - С. 52 с.

12. Лобатый, А. А. Структурно-параметрическая нечеткая коррекция алгоритма фильтрации / А. А. Лобатый, А. С. Бенкафо, А. С. Абуфанас. - М.: Системный анализ и прикладная информатика, 2014. - С. 4-8.

References

1. Krasilchikov, M. N. Modern information technologies in problems of navigation and guidance maneuverable unmanned aerial vehicles / M. N. Krasilchikov, G. G. Serebryakov. - M.: FIZMATLIT, 2009. - 558 p.

2. Aleshin, B. S. Orientation and navigation of mobile objects: modern information technology / B. S. Aleshin, K. K. Ve-remeyenko, A. I. Chernomorsci. - M.: FIZMATLIT, 2006. - 424 p.

3. Balakrishnan, A. V. Kalman filtering theory: Per. from English. / A. V. Balaknishnan. - M.: Mir, 1988. - 1000 p.

4. Sinitsyn, 1 N. Kalman filters and Pugachev / I. N. Sinitsyn. - M.: University Book, 2006. - 640 p.

5. Pupkov, K. P. Methods of classical and modern automatic control theory. 5 t / K. P. Pupkov, N. D. Egupov. - M.: Publisher MSTU. NE Bauman, 2004. - 656 p.

6. Tikhonov, V. I. Markov Processes / V. I. Tikhonov, I. N. Mironov. - M.: Soviet Radio, 1977. - 450 p.

7. Pugachev, V. S. Theory of stochastic systems / V. S. Pugachev, I. N. Sinitsyn. - M.: Logos, 2004. - 1000 p.

8. Wentzel, E. S. Probability theory / E. S. Wentzel. - M.: Higher School, 1999. - 576 p.

9. Kazakov, I. E. Methods of optimization of stochastic systems / I. E. Kazakov, D. I Gladkov. - M.: Nauka, 1993. - 270 p.

10. Lobaty, A. A. Features of the Kalman-Bucy filter complexes in orientation and navigation / A. A. Lobaty, A. S. Ben-kafo. - Minsk: Reports BSUIR, 2013. - P. 67-71.

11. Lobaty, A. A. Evaluation of navigation parameters of the movable object in a multimode / A. A. Lobaty, A. S. Benkafo. -Minsk: Reports BSUIR, 2014. - P. 52-58.

12. Lobaty, A. A. Structural-parametric correction fuzzy filtering algorithm / A. A. Lobaty, A. S. Benkafo, A. S. Abufanas. -Minsk: System Analysis and Applied Computer Science, 2014. - P. 4-8.

Поступила 01.03.2016

LOBATY A. A., YACINA Y. F., AREFYEU N. N.

OPTIMAL ESTIMATION OF RANDOM PROCESSES ON THE CRITERION OF MAXIMUM A POSTERIORI PROBABILITY

Belarusian National Technical University

The problem of obtaining the equations for the a posteriori probability density of a stochastic Markov process with a linear measurement model. Unlike common approaches based on consideration as a criterion for optimization of the minimum mean square error of estimation, in this case, the optimization criterion is considered the maximum a posteriori probability density of the process being evaluated.

The a priori probability density estimated Gaussian process originally considered a differentiable function that allows us to expand it in a Taylor series without use of intermediate transformations characteristic functions and harmonic decomposition. For small time intervals the probability density measurement error vector, by definition, as given by a Gaussian with zero expectation. This makes it possible to obtain a mathematical expression for the residual function, which characterizes the deviation of the actual measurement process from its mathematical model.

To determine the optimal a posteriori estimation of the state vector is given by the assumption that this estimate is consistent with its expectation - the maximum a posteriori probability density. This makes it possible on the basis of Bayes 'formula for the a priori and a posteriori probability density of an equation Stratonovich-Kushner.

Using equation Stratonovich-Kushner in different types and values of the vector of drift and diffusion matrix of a Markov stochastic process can solve a variety offiltration tasks, identify, smoothing and system status forecast for continuous and for discrete systems. Discrete continuous implementation of the developed algorithms posteriori assessment provides a specific, discrete algorithms for the implementation of the on-board computer, a mobile robot system.

Keywords: mathematical model, state space, stochastic equations, optimization criteria, the expectation, the demolition of the vector, diffusion matrix.

Lobaty A. A.

Doctor of science, professor. In 2000 he established chair «Information systems and technologies» in Belorussian National Technical University, department of «International Institute of Distance Education». His research interests include algorithms, concepts, and architecture for digital signal processing systems. He has extensive consulting experience in control of unmanned aerial vehicles. He is author and coauthor of many papers in scientific magazines, conference proceedings, and a number of books. He has number of university and state awards for achievements in teaching and research.

Yacina Y. F.

Director of the State Research and Production Enterprise unmanned multipurpose complexes. A specialist in the field of research and development of control systems of mobile robotic systems for various application. Head of research and development project on the creation of unmanned aerial vehicles.

Nikolay Arefiev received the graduate degree in software engineering from the Belarusian National Technical University in 2014 and the Master's degree in system analysis and control of information processing in 2015. He is currently working on PhD degree program. His current research interests include control of unmanned aerial objects.