Эффективный выбор плотности распределения случайной сетки при решении задачи о слежении по азимуту методом Монте-Карло

УДК 519.6

Д.Г.Арсеньев, В.М. Иванов, Н.А. Берковский

эффективным выбор плотности распределения случайной сетки при решении задачи о СЛЕЖЕНИИ ПО АЗИМУТУ методом МоНТЕ-КАРЛо

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

Последнее время этой задаче посвящено значительное число работ [5, 6, 8], в которых она, как правило, решается при помощи частичных фильтров (particle filters), в частности, различными модификациями схемы bootstrap filter [5]. Отметим, что задача является нелинейной задачей оптимального оценивания. [3, 4].

В настоящей работе применяется непосредственное интегрирование и предлагается новый, более эффективный метод выбора плотности распределения точек при вычислении интегралов в соответствии со схемой выборки по важности (important sampling) [5] в байесовском оценивании. С помощью этого метода при достаточно реалистичных предположениях об измерительном шуме и при условии малости системного шума задача решается значительно быстрее и точнее, чем немодифицированными последовательными методами Монте-Карло или методом выборки по важности при обычном для байесовского оценивания выборе плотности распределения узлов сетки.

Постановка задачи

Рассматривается задача оптимального оценивания с дискретным временем. Пусть наблюдаемый объект с координатами (x, y) движется по плоскости в соответствии со стандартной моделью [6] вида:

x, = Фх k- +rw k, k e N,

Ф

( 1 t 0 0 '0,5t 0

0 1 0 0 1 0

, г=

0 0 1 t 0 0,5t

v 0 0 0 1 V 0 1

t - шаг по

времени, равный интервалу между измерениями.

Системный шум (м>х, wy) будем считать гаус-совским белым шумом с средним значением (0,0) и матрицей ковариации qE, где q > 0, Е —двумерная единичная матрица. Пусть наблюдатель находится в начале координат и совершает измерения хк, которые зашумлены и определяются формулой:

-1 Г Уку

Z = tan

+ vk,k eN,

где измерительный шум V — также гауссовский белый шум с нулевым средним и среднеквадра-тическим отклонением г. Кроме того, известно априорное распределение вектора состояния в «нулевой» момент времени, в данной статье это нормальное распределение со средним и ковариационной матрицей стЕ для начального положения и равномерное распределение в квадрате 5" = [ах,Ъх] х [ау,Ъу ] е Я2 для начальных компонент скорости. При этом считаем, что начальные координаты и скорости взаимно независимы.

Требуется, располагая рядом измерений гх.к = г2, ..., гк} в к-й момент времени и знанием априорного распределения вектора сот

стояния хо = (х, х, у, у)о , приближенно найти апостериорное среднее значение Ер( | ) [х0Л:]

случайного процесса

=

:\T 0:к ,

где

где xk = ( - x y, y )1,

w k = (wx , Wy )T ,

Х0:к = X1, ..., Хк }, а Р(Х0:к I 21:к ) - апостериорнЗЯ

плотность. Здесь придерживаемся обозначений, принятых в теории байесовского оценивания [5].

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

x

Рис. 1. Постановка задачи слежения по азимуту

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

Отметим, что результатом байесовского оценивания является апостериорная средняя траектория лодки, вычисленная по итогам полученных измерений (2) с учетом соотношений (1). Эта траектория не обязана совпадать с истинной траекторией лодки, но при достаточном числе измерений и небольших шумах весьма близка к ней.

В момент времени к имеем байесовскую оценку для средней траектории в виде [3, 5]:

без системного шума с известной начальной точкой, т. е. модель вида

хк = фхк-1 , х0 =

где (х0, у0) - фиксированная начальная точка, а вектор скоростей (х, у) равномерно распределен в квадрате S. В этом случае интегралы в (3) будут двукратными, а знаменатель будет равен:

j P(zi:k 1 X0:k) P(X0±)dX0± =

= j|Sf1 П(2nr2)-0'5 expx

E

P(x0:k |zl:k H 0:k

K:k ] =

j X0:kP(Z1:k 1 X0:k )P(x0:k )dx0:k (3) j P(zi:k 1 X0:k )P(x0:k )dx0:k

i=1 0

z, - tan

v

y + yit

\\

^ X Xit j j

(2r2)-

dXdy.

где р(х0.к) - априорное распределение траекторий, определяемое начальным распределением вектора состояния и распределениями системных шумов. Интегралы в общем случае имеют кратность к + 4 и при большом количестве измерений к их приближенное вычисление является проблемой.

Выбор плотности распределения для скоростей

Вначале исследуется вспомогательная задача

График подынтегральной функции из (5) показан на рис. 2

Как видно из рис. 2, эта функция имеет единственный резкий максимум. Хотя ее график и напоминает внешне нормальное распределение, в общем случае она асимметрична и имеет ненулевой эксцесс.

Предположим, что мы приближенно вычисляем числитель и знаменатель в (3) методом Монте-Карло и используем дельта-метод [2, 7] для оценки частного. Действуя классическим методом important sampling, для генерации случайной сетки следовало бы использовать

х 10 '

7 -10

Рис. 2. Типичный график апостериорной плотности для скоростей

равномерное распределение в квадрате. Однако численные эксперименты и теоретические соображения показывают, что для уменьшения количества итераций очень хорошим выбором является распределение, пропорциональное всей подынтегральной функции в (5). К сожалению, оно, по всей видимости, недоступно для эффективного моделирования, например, в методе отбора затруднительно найти более или менее точную мажоранту, и подавляющее большинство генерируемых точек не проходит отбор, соответственно, время счета сильно увеличивается. Классический же метод important sampling

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

Решаем систему неравенств относительно координат точки (х, у) е £ :

(2пг2) 0 5 exp - z,. - tan-

К x, y )е 5

(

0

y + yit

x0 + Xit

(2r2)-

yj

> c, i = 1, 2...k

где с > 0 - некоторое достаточно малое число (в численных примерах в настоящей работе оно принято равным 0,05). Имеем всего к + 2 неравенств. Если предположить, что

х >-х\иу\ г = 1,2... к ,

то из (6) можно прийти к системе неравенств

ах а . „ . Ь.х а

ах + - У° < у < Ьх + - у0 г = 1,2... к п и ,

(х, У) е £

где аг = гап( - г у/ 21п(с(2пг2)0,5)-1) и

Ь, = гап( г, + Гу/ 21п(с(2пг 2)а5)-1).

Ввиду линейности относительно х и у двойных неравенств системы (8), множество решений системы (8), если оно существует, будет представлять собой некоторый многоугольник. На рис. 3 изображена типичная ситуация в случае когда все положения объекта наблюдатель видит под углом от нуля до 90 градусов. При этом скорости равномерно распределены в квадрате £ = [-10,10] х [-10,10], число измерений равно 15, из-за этого так много прямых на рисунке. Для наглядности показана только часть квадрата £, в которой лежит область решений системы неравенств (8).

Рис. 3. Многоугольная область решений системы неравенств

Предположим, что нам известны координаты достаточного числа точек, удовлетворяющих системе неравенств (8) (технически они легко находятся с помощью простой и быстрой подпрограммы). Располагая этим «облаком точек», как принято говорить в статистике, мы можем найти выборочное среднее х = (х, у), а также ковариационную матрицу £ х этого облака.

Далее, в качестве плотности для генерации случайной сетки при вычислении методом Монте-Карло интегралов из (3) возьмем усеченное нормальное распределение вида:

Рх

где х = (х, у), X, - характеристическая функция квадрата, к > 0 - параметр, а ^х, х, ) здесь и далее обозначает функцию плотности нормального распределения с аргументом х, средним х и ковариационной матрицей . В отличие от распределения, пропорционального выражению

П exp

с

f о

z, - tan-

У + yit

w

V

y X + Xlt y^

2\ -0,5

из (5),

это распределение моделируется эффективно. В то же время оно сгущает случайную сетку в районе точки максимума апостериорной плотности (рис. 2 и 4).

Управляя параметрами с и к, можно еще больше уменьшить дисперсию. Отметим, что в случае,

когда многогранник области решений системы (8) имеет вид треугольника с прямым углом в вершине квадрата, (что достаточно редко), то вместо нормального распределения лучше использовать асимметричные параметрические распределения, поддающиеся эффективному моделированию. В данной статье эти ситуации не рассматриваются.

Теперь перейдем к решению задачи (3) при условии отсутствия системных шумов. В этом случае интегралы в (3) становятся четырехкратными, а знаменатель имеет вид:

j P(zi:k \ X0:k )P(X0:k )dX0:k = = j| S\П (2nr2)-°,5expx

f

z - tan-

У + yit

X + Xlt

(2r2)-

x N(x, x0, Ex )dXdydx,

где x = (X, y), x0 = (X0, y0 ) и px (x) - априорная плотность распределения вектора координат объекта. Принимая во внимание изложенные выше соображения, в качестве плотности для генерации случайных точек рационально взять плотность

Px

где Px (x) определена равенством (9). Заметим, что в случае классической схемы important

i=i

0,7^ 0,6^

0,5 -,

Рис. 4. График подобранного нормального распределения при к = 1, соответствующий плотности из рис. 2

sampling [5] соответствующая плотность имеет вид

,1-1

p (x) = N( x, x0, S x )| S\

Далее, экстраполируем эти же соображения на общую задачу с системным шумом.

Отметим, что для того чтобы оценить траекторию объекта, достаточно оценить начальное положение, начальную скорость, а также процесс системного шума ^^к - 2, ..., к}. По оцененным значениям системных шумов, легко восстановить траекторию, пользуясь соотношениями модели (1). Введем вектор X ^ = (х0, х0, у0, у0, w2, ..., wk), подразумевая, что каждое вхождение w. обозначает две компоненты соответствующего системного шума, т. е. вектор xW имеет размерность 2к + 4 . Дополнительно введем вектор 1, = (г - 0,5, г -1,5,..., 0,5)г и обозначим как w'х и w'y векторы, состоящие из первых и вторых компонент векторов w., где у меняется от единицы до г, а г - от единицы до к. Требуется вычислить байесовскую оценку вектора xW по формуле:

>[Х-]- Гр^кКх ,

причем

P

z - tan "

y + yit + tw y l , x + xit + twi l.

•

2\ -0.5

Априорная плотность p(xW ) является произведением равномерной плотности для начальных скоростей, нормальных плотностей для начальных координат и 2k шумов. Таким образом, требуется вычислить 2k + 4 интеграла размерности 2k + 4. Действуя по классической схеме important sampling, в качестве плотности для генерации случайной сетки следует взять функцию

px0 « ) = N(x, x0, S x ) px (х)П N(w i ,(0,0), qE). (15) w i=1 Вычисляя интеграл предлагаемым методом, для генерации случайной сетки берем плотность

px0 (xW ) = N(x, x0, S x )| S\П N(w i ,(0,0), qE). (16) w i=1 Таким образом, (15) и (16) различаются только способом генерации случайных сеток для начальных скоростей, однако, как показано ниже в разделе численных экспериментов, при малых системных шумах это отличие имеет решающее значение. Отметим, что перспективной представляется возможность адаптивного управления параметрами c и h в предлагаемой схеме с привлечением теории, развитой в [1].

Результаты численного эксперимента

Ниже приведены результаты решения двух

V

i—1

Таблица 1

Сравнение методов при отсутствии системного шума

Метод Предлагаемый метод с эффективно выбранной плотностью Метод выборки по важности Частичный фильтр

Число итераций 10 000 380 000

Время, с 4 95 368

задач. Первая задача - это задача о слежении по азимуту при условии нулевого системного шума. Среднеквадратические отклонения модели равны r = 1°, о = 10 м, шаг по времени t = 5 с, параметр h = 1. Количество измерений равно 15. Начальные скорости, измеряемые в м/с, равномерно распределены в квадрате S = [-10,10] х [-10,10] е R2.

Байесовская оценка моделируемых траекторий вычислялась тремя методами: классическим методом выборки по важности для байесовского оценивания с плотностью случайной сетки, определяемой формулой (12), предлагаемым методом с эффективно выбранной плотностью, определенной по (11), и при помощи частичного фильтра (немодифицированная схема bootstrap filter [5]), использующего 20000 частиц.

Расчет проводился на доверительном уровне 95 % с точностью 1 м для координат начального положения и 0,1 м/с для начальных скоростей.

Результаты моделирования показаны на рис. 5 и в табл. 1, из которых видим, что средние траектории, полученные предложенным методом и методом important sampling, совпадают, как и следовало ожидать. Кроме того, они ближе к истинной траектории, чем траектория, оцененная методом bootstrap filter. Поскольку схема bootstrap filter использует дискретные аппроксимации априорной и апостериорной плотности на каждом шаге по времени, ошибка накапливается. Причина столь долгого выполнения bootstrap filter - необходимость моделирования дискретного распределения с большим числом частиц на каждом шаге при

300 400 500 600 700 800 900 1000

Рис. 5. Задача без системного шума

1 - траектория объекта; 2 - оценка методом выборки по важности; 3 - оценка предложенным методом; 4 - оценка, использующая частичный фильтр

Сравнение методов при наличии малого системного шума

Таблица 2

Метод Предлагаемый метод с нулевым системным шумом Предлагаемый метод с учетом системного шума Метод выборки по важности Частичный фильтр

Число итераций 10 000 620 000 6 770 000

Время, с 6 758 6 970 5 601

«перевыборе» частиц. Именно на эту процедуру программа тратит основное время. Из табл. 1 видим, что предложенный метод значительнее эффективнее двух других как в отношении времени работы, так и в отношении количества итераций при сравнении с important sampling.

Вторая задача - это задача о слежении по азимуту с малым системным шумом. Значения параметров те же, что и в предыдущей задаче, но h = 8, это оказался наиболее оптимальный вариант. Дополнительно вводился системный шум со среднеквадратическим отклонением q = 0,1 м/с. Расчет проводился на доверительном уровне 95 % с точностью 2 м для координат начального положения, 0,1 м/с для начальных скоростей и 0,01 м/с для системных шумов на каждом шаге.

Для оценки точности частного интегралов использовался дельта-метод [7].

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

Результаты моделирования показаны на рис. 6 и в табл. 2, из которых видно, что траектории, оцененные предлагаемым методом и методом выборки по важности, снова совпадают. Менее точно повторяет истинную траекторию оценка, получен-

900

800

700

600

500

400

300

200

100

2,3 5 \ ..................... ..................... ! I

1 ^О^

; Начальная ^^ | точка| i i

400

500

600

700

800

900

1000

Рис. 6. Задача с малым системным шумом 1 - траектория объекта; 2 - оценка методом выборки по важности; 3 - оценка предложенным методом; 4 - оценка, использующая частичный фильтр; 5 - решение предложенным методом с нулевым

системным шумом

ная предлагаемым методом без учета системного шума, и наиболее неточны результаты, полученные с помощью частичного фильтра.

Однако отметим, что время работы как предлагаемого метода, так и метода выборки по важности резко увеличивается, т. к. в этом случае мы переходим от четырехкратных интегралов к 34-кратным. Можно выдвинуть тезис о целесообразности пренебрежения малыми системными шумами, т. к. выигрыш по времени в этом случае будет очень большой (см. табл. 2), а проигрыш в точности невелик (рис. 6). Наконец, отметим, что при больших системных шумах (составляющих порядка 10 % от характерной скорости) преимущество предлагаемого метода над классическим методом important sampling теряется, а замена решения задачи с системным шумом на решение задачи при его отсутствии приводит к существенным неточностям.

В статье предложен новый, эффективный метод выбора плотности распределения узлов случайной сетки для вычисления байесовской оценки для задачи о слежении по азимуту. Он основан на рациональном выборе плотности распределения точек для генерации случайной сетки методом Монте-Карло. Выбирается плотность, которая более близка к оптимальной плотности, чем априорная плотность распределения вектора состояния, и в то же время доступна для эффективного моделирования.

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

список литературы

1. Арсеньев, Д.Г. Адаптивное управление в стохастических методах вычислительной математики и механики [Текст] / Д.Г. Арсеньев, В.М. Иванов, М.Л. Кореневский; 2-е изд., испр. и доп. - СПб.:Наука, 2008. -423 с.

2. Арсеньев, Д.Г. Навигация по расстояниям до точечных ориентиров адаптивным методом существенной выборки. [Текст] / Д.Г. Арсеньев, В.М. Иванов, Н.А. Берковский // Научно-технические ведомости СПбГПУ -2011. -№ 1 (115). -С. 81-86.

3. Cтепанов, О.А. Введение в теорию оценивания. [Текст] / О.А. Степанов // Основы теории оценивания с приложениями к задачам обработки навигационной информации: Ч. 1. -СПб.: ГнЦ РФ ЦНИИ «Электроприбор», 2009. -496 с.

4. Stepanov, O.A. Investigation of Linear Optimal Estimator [Text] / O.A. Stepanov, A.B. Toropov // Proc. of

17-th World Congress. -Seoul, July 6-11. -P. 2750-2755.

5. Doucet, A. Sequentiual Monte-Carlo methods in practice [Text] / A. Doucet, N. Freitas, N. Gordon. - N.Y.: Shpringer-Verlag, 2001. -581 p.

6. Gordon, N.J. Novel approach to nonlinear/non-Gaussian Bayesian state estimation [Text] / N.J.Gordon, D.J. Salmond // IEE Proc. -Apr. 1993. -Vol. 140. -№ 2. -P. 107-113.

7. Geweke, J. Bayesian inference in econometric models using Monte-Carlo integration [Text] / J. Geweke // Econometrica. -Nov. 1989. -Vol. 57. -№ 6. -P. 1317-1339.

8. Kostas E. Bekris. Evaluation of Algorithms for Bearing-Only SLAM// [Электронный ресурс] / Kostas E. Bekris, Max Glick, L.E. Kavraki. -2006. -Режим доступа: http://www.cse.unr.edu/~bekris/papers/bearing_only_ slam.pdf

УДК 681.3.06

Б.Г. Ильясов, И.В. Дегтярева, Е.А. Макарова, Т.А. Карташева

интеллектуальные алгоритмы принятия решений при управлении инвестиционным процессом макроэкономической системы

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