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

УДК 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)

Ф

( 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,

(2)

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

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

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

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

= (x, x, y, y)

:\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, Ъ У \ У), (4)

где (х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

(5)

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

(6)

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

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

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

ах а . „ . Ь.х а

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

(х, У) е £

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

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

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

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

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

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

Рх(х) = ^х,х,к•£х)-X, , (9)

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

П exp

с

f о

z, - tan-

У + yit

w

V

y X + Xlt y^

(2r2)

2\ -0,5

из (5),

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

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

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

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

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

(10)

f

z - tan-

У + yit

X + Xlt

(2r2)-

x N(x, x0, Ex )dXdydx,

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

Px (x) = N(x, x0, Ex )A (x), (11)

где 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\

(12)

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

Отметим, что для того чтобы оценить траекторию объекта, достаточно оценить начальное положение, начальную скорость, а также процесс системного шума ^^к - 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 по формуле:

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

причем

P(z1:к | xW) = ТГ(2nr2)-0'5 expx (14)

z - tan "

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

•(2r2)

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

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

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

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