Применение весового метода МонтеКарло к исследованию редких событий в динамических системах, возбуждаемых белым шумом Текст научной статьи по специальности «Математика»

УЧЕНЫЕ ЗАПИСКИ ЦАГИ Том XXVII 1996

№3-4

УДК 629.735.33.015.073

ПРИМЕНЕНИЕ ВЕСОВОГО МЕТОДА МОНТЕ-КАРЛО К ИССЛЕДОВАНИЮ РЕДКИХ СОБЫТИЙ В ДИНАМИЧЕСКИХ СИСТЕМАХ, ВОЗБУЖДАЕМЫХ БЕЛЫМ ШУМОМ

А. А. Дынников

Предлагается метод исследования статистических характеристик динамических систем, возбуждаемых белым шумом, являющийся модификацией весового метода Монте-Карло. Данный метод ориентирован на решение практических задач управления полетом летательных аппаратов, требующих • ' _й оценки вероятностей порядка 10 + 10 (например, исследование систем

автоматической посадки). Приведены формулы для определения весовых коэффициентов для различных вариантов преобразования исходных возмущений. Рассмотрены примеры применения предлагаемого метода и показано, что потребное число реализаций при оценке вероятностей порядка

10~7 + 10-11 может быть сокращено в 10^ + 105 раз.

1. Постановка задачи. Рассматриваемся динамическая система, описываемая дифференциальными уравнениями

х(і)'єЯп, ?(/) єЛк, / є[0,Г], Т = сопяг, х(0) = х0,

(1)

где х (0 — вектор фазовых координат динамической системы, £(*) —

вектор некоррелированных гауссовских белых шумов с единичной интенсивностью, g — матрица соответствующей размерности. Стохастическая система (1) понимается в смысле Стратоновича [1]. Подобные уравнения часто используются для решения задач динамики полета, например, при определении вероятности летного происшествия при автоматической посадке летательного аппарата. Согласно существующим нормативам, эта вероятность не должна превышать КГ6 * 10-8.

В последнее время был разработан ряд методов [2], [3], позволяющих находить асимптотические оценки вероятности для подобных задач. Реализация этих методов встречает трудности, связанные с необходимостью решения сложной задачи оптимального управления для систем большой размерности. Метод Монте-Карло в этом отношении является более простым, однако к таким задачам практически не применяется из-за большого потребного числа реализаций. Известно, что использование весовых методов Монте-Карло позволяет уменьшить потребное количество реализаций и снизить дисперсию получаемой оценки. Эти методы широко используются для решения многих задач, связанных с вычислением статистических характеристик различных систем [4]. Ниже предлагается метод, являющийся модификацией весового метода Монте-Карло для задач, включающих в качестве возмущений белый шум.

Рассматривается задача о нахождении среднего значения П функции П(х), определяемой соотношением

П°(х), Ф(х)>Ф,,

^ ’ 0, Ф(х)<Ф*.

Предполагается, что вероятность выполнения неравенства ф(х(7)) £ Ф* мала, т. е. функция П(х) редко принимает ненулевые значения, Частным случаем этой задачи является задача об определении вероятности Р = р(ф^х(Г)) > Ф„) превышения заданной функцией

ф(х(Т)} фазовых переменных некоторого заданного уровня Ф*:

ф(х(7’)) > Ф*. (2)

Эта задача эквивалентна задаче о нахождении среднего значения при а0(х) = 1. Ниже рассматривается задача о нахождений среднего значения П.

Предполагается, что доля траекторий, для которых выполнено неравенство (2), мала. Это возможно в двух случаях: если область в Яп , в которой выполняется неравенство (2), имеет малый объем, либо эта область велика, но расположена вдали от конца невозмущенной траектории, из-за чего плотность вероятности в этой области мала и интеграл от нее по всей области существенно меньше единицы. Будем рассматривать только второй случай, предполагая, что увеличением дисперсии х(Г) можно добиться значительного увеличения Доли событий,

для которых неравенство (2) выполнено.

2. Метод определения среднего значения. Рассмотрим следующую схему определения среднего значения (модификация весового метода Монте-Карло):

а) выполняется N однотипных реализаций, на каждой из которых интегрируется система (1), причем на вход динамической системы

на 1-й реализации вместо /-Й реализации вектора белого шума §*(/) подается величина £*(?) (здесь и далее верхний индекс означает номер реализации), определяемая соотношением

т

■ "■ -0 ' .. -где К (і, в}— кусочно-непрерывная матричная функция, заданная на

квадрате [О,Г] х [О,Г]. Эту функцию будем называть преобразующей

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

В результате интегрирования системы определяется соответствующее значение Ог функции О.;

б) полученные значения П' суммируются для определения оценки искомой величины П с весовыми коэффициентами р':

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

Покажем, что для этого весовые коэффициенты Должны определяться соотношением

где С — нормировочный коэффициент, определяемый из условия равенства единице среднего значения весовых коэффициентов р = 1, штрих (') обозначает транспонирование. Вопрос о вычислении нормировочного коэффициента б будет рассмотрен ниже. Равенство р = 1

Для доказательства соотношения (5) рассмотрим конечно-мерный случай, в котором белый шум представлен в виде ступенчатой функции

(4)

V |_о

Л , (5)

/

является необходимым условием несмещенности оценки (4). Это условие легко получить, представив, что П = О.1 & 1.

іде х = Т/М— ширина одной ступеньки, М — число ступенек, а величины — независимые случайные векторы, состоящие каждый из к независимых случайных величин, распределенных по нормальному закону С нулевым средним и дисперсией 1/т. Вектор I = (іі распределен по нормальному закону с плотностью распределения

=......

. (2я) 2 № 2

(6)

где Щ — корреляционная матрица вектора \ размера кМ х кМ, равная .

(I л

- О - О т

о - ••• О

х

О о - -

/ -г+- единичная матрица размера к х к.

_ /

Компоненты вектора С = } » являющегося ступен-

чатым приближением Г\м{*) функции Г|(/) , определяются соотношением

Су = ^ у + Л у, =

т=1

где К]т — элемент матрицы АГ = {^ут} , определяемый по формуле

1 п і

и-ъ

1Г$ V

I ^(а,р)ф

г/а.

Отсюда следует, что вектор С, распределен по нормальному закону с плотностью распределения

где корреляционная матрица вектора С определяется соотношением

^ =Л^+± (ЯК'+*'+*).

Условие несмещенности оценки О может быть получено следующим образом: с одной стороны, £} = с другой —

после замены | на £ и введения весовых коэффициентов 0=1 ЦС)^г(С)Рь(^)г^ * Заменяя в первом из этих интегралов

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

Подставляя в это равенство плотности вероятности (6) и (7), получаем, что весовой коэффициент должен равняться

= бехр(-|[с{^-1-^-1}с]) =

= £ехр+ к.+ *'*}!]) = = о ехр^- ^ + 2х\'\|],

где 0 = 2 — величина, не зависящая от £ .

Возвращаясь к функциям тц^(/) и £,м (*), получаем

р, = (?ехр

Устремляя М к бесконечности, пйлучаем соотношение (5).

В заключение отметим, что в тех случаях, когда определение нормировочного множителя <2 затруднительно, в (5) можно использовать

~ 1 * •

его оценку 0 = —У*, р1 . При этом выражение (4) принимает вид /=1

3. Способы выбора преобразующей формы. Известно, что дисперсия оценки при использовании весового метода Монте-Карло сильно зависит от выбора преобразования для плотности вероятности [4]. В рассматриваемом методе она определяется способом выбора преобразующей формы (3). В общем случае преобразующую форму следует выбирать из условия минимума дисперсии получаемой оценки

Da = (п* р1' - n)2 -» inf, (8)

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

CljN

уменьшая Dq и увеличивая N. Последнее, однако, не всегда возможно. Например, в задачах оценки малых вероятностей, потребное число реализаций N может оказаться слишком большим, и следует попытаться уменьшить величину Dq. Следует отметить, что минимум в (8)

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

4. Случай линейной системы с постоянными коэффициентами и фиксированным временем. Приводимый ниже пример иллюстрирует возможность выбора преобразующей формы в соответствии с условием (8).

Рассматривается динамическая система, описываемая линейным дифференциальным уравнением:

x(t) = A{t)x(t) +

х(0) = 0, x(t) є Rn, §(it) eRk, t є[0,Г], T = const,

где A(t) и g(t) — матрицы соответствующей размерности.

«йг»=

і, Ф{х(Т))>Ф,> ^x) = h,x

О, Ф(х(Т))<Ф., -Ч' •

В интегральной форме уравнение (9) имеет вид

'■ > :' ■

х(і) = | ИГ(/,т)я(>)£(т)Ж,

О

где И^,т) — решение матричного дифференциального уравнения

dt

Отсюда следует, что

где q{i) — импульсная переходная функция:

q(t) = h'W(T,t)g(t).

(П)

(12)

Из соотношения (11) следует, что Ф распределено по нормальному закону с нулевым средним и дисперсией |Ы|2:

Рф(х) =

ÆHI

ехр

m

Здесь и далее под нормой ||#| функции д понимается норма

(т

и-

vo >

Остановимся на вопросе выбора преобразующей формы. При любой реализации £(/) добавок ri(/) в (3) может быть представлен в виде

1 Т

П(0 = “?(*) + Р(0 » гДе а = ~2 f ^'{t)y\{t)dt, р(/) = л(') ~ «q(t). Из (11)

. , N1 0

следует, что слагаемое p(i) не будет влиять на распределение Ф при

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

няться нулю. Немного более сложными рассуждениями, которые здесь не приводятся, можно показать, что матричная функция АГ(/,б), удовлетворяющая условию (8), с точностью до множителя о должна иметь

ввд

м

Исследуем зависимость дисперсии оценки от параметра а.

Математическое ожидание оценки П, представляющее собой вероятность события Ф > Ф*, очевидно, не зависит от а и равно

00 X

О

1 е.-Т

/

где Л = -тру. После замены £(/) на £(/) в соответствии с (3) имеем

Ф = |= (1 + ст)|, о о

причем плотность распределения величины Ф теперь определяется соотношением

?ф(*)

1

Ґ

ехр

I 2|ИГ(1^)

Дисперсия оценки

ЛП-

п2

= і [од

■со .

(1 + о)3

£ф(х)

гг

Рф(х)сЬс и

2 ст + 4ст +1

Ґ -у Л

Я2

12(1 + а)2)

-1

(13)

Д£) (Л.с)

График зависимости функции 5\р,Щ = —Щд)' 07 ст ПРИ Различ~

ных значениях Л приведен на рис. 1. Эта функция показывает, насколько может быть уменьшено число реализаций, потребное для оценки П, так как это число пропорционально Я(р, Л):

N

,, ¿М)

Рис. 1. Зависимость функции £(а, Л) от о при различных значениях Я

Для обычного метода Монте-Карло ст = 0 , N = ==■. При ст = О

функция *У(ст,7?) равна единице и убывает с ростом ст . При ст -» со она

нео1раниченно возрастает. Выбирая ст так, чтобы *У(ст,Л) было

минимально, можно значительно сократить потребное количество реализаций.

Так, для П » 2,8 • 10-7 (Л = 5) возможно, выбирая ст и 4, уменьшить потребное число реализаций примерно в Ю5 раз. Как видно из рис. 1, функция ^(ст,Л) в окрестности минимума меняется слабо,

поэтому точное его определение не обязательно.

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

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

Дисперсия оценки (8) зависит, с одной стороны, от разброса значений оцениваемой величины, с другой — от дисперсии весовых коэффициентов. Поскольку доля реализаций, для которых Л1 * О, невелика, уменьшить дисперсию оценки (8) можно, увеличивая долю таких реализаций. Отсюда следует, что преобразующую форму нужно выбрать так, чтобы увеличить отклонения фазовых переменных системы (1). Этого можно достичь, изменяя спектральную плотность исходных возмущений.

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

Целенаправленная деформация спектральной плотности может осуществляться, например, формирующим фильтром

Л'(0) = Т10,

что соответствует частному случаю преобразования вида (3) !

t

У'(/) = |*(/,е)^'(е)й/е, о

где ■£■(*, 0) удовлетворяет уравнениям

к (в,в) = <?(е).

Отметим также, что в данном случае преобразование возмущений осуществляется не на квадрате [0,7] х [0,Г], а на треугольнике

0 £ 0 £, t й Т, что более удобно для практической реализации статистического моделирования, так как для преобразования возмущений не требуется знать «будущих» (0 > /) возмущений.

В качестве примера применения такого подхода рассмотрим самый простой случай уравнения первого порядка. Требуется опреде-

(14)

лить вероятность превышения величиной х(Г) уровня х*. Рассмотрим дифференциальное уравнение

х(/) = -х(/) +

і є[-оо,Г],

(15)

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

В соответствии С приведенными выше рассуждениями ДЛЯ Г|(/) был взят формирующий фильтр вида

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

Были выполнены серии испытаний по 40 ООО реализаций для различных значений параметра сі при Т = 10 (выбор именно такого значения Т объясняется тем, что для практических задач динамики полета характерное отношение времени Т в (1) к времени автокорреляции равно ~ 10 т 30). Полученные результаты представлены на рис. 2—5.

Щ = -л(0 + <%=(?).

/ > 0, Л(0) = 0.

вценибаемая вероятность

е

□

° в в О

□

10

□

о

□

Рис. 2

о — весовой метод Монте-Карло, 40 000 реализаций, й * 2,0; • — вероятность реализаций в процессе моделирования без учета весовых коэффициентов; О — истинное значение вероятности

»

да'

■* иґ 10* 10* 10~Я

------------—і І I ---

"—1 'Т І~ I і" I..................."I .....Г' “Т" і І Iі

....................................................

...........................

вв»п_

0В,

о

оно

о

В

в

Рис. 3 Обозначения такие же, как на рис. 2 ОчеШимая кратность

- Рис. 4. Оценка величины относительной погрешности

для различных значений параметра й

На рис. 2 и 3 представлены оценки вероятностей, полученные при значениях ё = 1,5 и 2,0 соответственно. По оси X отложено истинное значение вероятности, вычисляемое как

^)-^гН-тГ

По оси У отложена оценка вероятности, полученная с помощью предлагаемого метода при соответствующем значении параметра <і, а также истинное значение вероятности и вероятность реализации события х(Г) а х* в процессе моделирования. Для значения <1 = 1,5 хорошее совпадение оценки с истинным значением наблюдается вплоть до значений Р * 10-9, а для г/ = 2,0 — вплоть до значений Р я ДО-11. На рис. 4 показана оценка величины относительной погрешности, приходящейся на одно испытание, для различных значений параметра сі:

р', х1 2> X,, 0, X* < X,.

Рис. 5. Оценка потребного числа реализаций для различных значений параметра й

Учитывая, что с ростом числа реализаций погрешность оценки вероятности убывает как ст « с/л/Ж, потребное число реализаций может быть оценено как N »ст2. Оценка этой величины для различных значений параметра й приведена на рис. 5. Из него следует, что даже для значений вероятности Р «10-11 потребное число реализаций составляет N «105, в то время как при применении обычного метода

Монте-Карло потребовалось бы порядка 1011 реализаций. Таким образом, для данного примера предлагаемый метод позволяет сократить число испытаний на пять порядков.

Работа выполнена при поддержке РФФИ.

ЛИТЕРАТУРА

1. Стратонович Р. Л. Новая форма записи стохастических интегралов и уравнений // Вестник МГУ. Серия мат. и мех.— 1964, № 1.

2. Кузьмин В. П., Ярошевский В. А. Асимптотические оценки вероятностей больших случайных отклонений фазовых координат динамических систем от средних значений // Ученые записки ЦАГИ.— 1990.

Т. XXI, № 3.

3. Дынников А. А. Асимптотическая оценка вероятностей больших случайных отклонений фазовых координат динамической системы при нали-

. чии фазовых ограничений // Ученые записки ЦАГИ.— 1994. Т. XXV, № 4.

4. М и х а й л о в Г. А. Оптимизация весовых методов Монте-Карло,—

М.: Наука.— 1987.

Рукопись поступила 27/1Х1994 г.