close

Вход

Забыли?

вход по аккаунту

?

Применение весового метода МонтеКарло к исследованию редких событий в динамических системах возбуждаемых белым шумом.

код для вставкиСкачать
УЧЕНЫЕ
ЗА.ПИСКИ ДА.rи
XXYII
Том
.NilЗ-4
1996
УДК 629.735.33Ш5.073
ПРИМЕНЕНИЕ ВЕСОВОГО МЕТОДА МОНТЕ-КАРЛО
К ИССЛЕДОВАНИЮ РЕДКИХ СОБЫТИЙ В
'ДИНАМИЧЕСКИХ СИСТЕМАХ, 'ВОЗБУЖдАЕМЫХ
БЕЛЫМ ШУМОМ
А. А. Дыннu1W8
Предлaraeтcи мeroд исследования: статистических характеристик дина­
мических сиcreм, возбуждаемых белым шумом, ивляющийся модификацией
весового :метода Монте-Карло. данный метод ориентирован на решение
прaJcrИЧCCXИX задач управлении полетом летательных аппаратов, требующих
,
'
оценки вероятнocreй порядка
'-б
+ 10
-8
(например, исследование сиcreм
автоматической посадки). Приведены формулы для определении весовых ко­
10
эффициентов для ра3ЛИЧНЫХ вариантов npeoбразования: исходных Возмуще­
НИЙ. Рассмотрены npимеры npименения: предлагаемого метода и показано,
что
потребное число реализаций при оценке
вероятноcreй порядка
10
1.
-7
-113
5
может быrь сокращено в 10 + 10
+ 10
Постановка задачи.
раз.
РассматрищеJCЯ динамическая с,истема,
описываемая дифференциальными уравнениями
X(t), = f{x(t),t) + g{x(t),t~(t),
x(t) 'Е я" , ~(t) Е Rk , t Е [О, Т],
х(о)
где х (1) -
Т = const,
(1)
= хо,
вектор фазовых координат динамической системы, ~(t)­
вектор некоррелированных гауссовских
интенсивностью,
стическая
}
система
белых шумов
с единичной
матрица, соотве;rcтвующей размерности. Стоха-
g -
понимается
(1)
в
смысле
Стратоновича
[1].
Подобные уравнения ,часто используются для решения задач динамики
полета,
например,
происшествия при
Согласно
.при
определении
автоматической
существующим
посадке
нормативам, ,эта
вероятности
летного
летательного
аппарата.
вероятность" не
должна
превышать 10-6 + 10-8 .
131
в последнее время был разработан ряд методов
[2], [3],
позво­
ляющих находить аСИМIП'OТИЧескИе оценки вероятности для подобных
задач. Реализация этих методов встречает трудности, связанные с не­
обходимостью решения сложной задачи оптимального управления для
систем большой размерности. Метод Монте-Карло в этом отношении
является более простым, однако к таким задачам практически не при­
меняется из-за большого потребного числа реализаций. Известно, что
использование весовых методов Монте-Карло позволяет уменьшить
потребное количество реализаций и снизить дисперсию получаемой
оценки. эти методы широко используются для решения многих задач,
связанных с вычИслением статистических характеристик,РЗ3.ЩlЧНЫХ си­
стем
[4].
Ниже npе.м:агается метод, являющийся МОДИфИIOlЩlейвесово­
го метОда . Монте-Карло для задач, включающиХ в :качестве вО~мущений
белый шум.
Рассматривается
задача
о нахо:ж;цении среднего
значения
n
функции п(х) , определяемой соотношением
n(Х) = Ino(x), Ф(Х) ~ ф.,
О,
ПреДполагается,
что
Ф(х) < Ф •.
вероятность
выполнения
неравенства
Ф(х(Т)) ~ Ф. мала,· т. е. функция n(х) редко принимает ненулевые
значения, Частным случаем этой задачи является зад~ча об определе­
нии веро~ности Р = Р(Ф(х(Т)) ~ Ф.) превышения заданной функцией
Ф(Х(Т)) фазовых переменных некотороro заданного уровня Ф.:
Ф(х(Т)) ~ Ф•.
(2)
эта задача эквивалентна задаче о нахо:ж;цении среднеro значения при
По(х);;
1.
Ниже рассматривается задача о нахо:ж;цении' среднего значе-
ния П.
Предполагается, что доля траекторий, для которых выполнено не-
равенство
(2),
мала. Это возможно в двух случаях: если область в
которой выполняется неравенство
(2),
Rn ,
в
имеет малый объем, либо эта
область велика, но расположена вдали от конца невозмущенной траек­
тории, из-за чего плотность вероятности в этой области мала и инте­
трал от нее по всей области существенно меньше единицы. Будем рас­
сматривать только второй случай, предполагая, что увеличением дис­
персии· Х(Т) можно добиться значительноro увеличения доли соБЫтий,
для которых неравенство
2.
(2)
выполнено.
Метод определения среднего значенни. РассмОТрим следующую
схему определения среднего значения (модификация весового метода
Монте- Карло):
а) выполняется
N
однотипных реализаций, на ка:ж;цой из кото­
рых интетрируется система
132
(1),
причем на вход динамической системы
на i-й реализации вместо i-й реализации веЮ'Ора белого шума
~i (/)
подается величина ~! (t) (здесь и далее верхний иНдекс означает номер
реализации), определяемая соотношением
т
f
;
~; (/) == 111 (t) + ~i (/),' 11 (t) = K(/,e~i (e~,
(3)
о
где
K(t,e)-
кусочно-непрерывна:я матричная ФУНКЦИЯ, заданная на
квадрате [О, Т] х [О, Т]. эту ФУНIЩИЮ будем называть преобразующей
формой. как будет по:казано ниже, выбор этой Функции определяет
эфФекгивность рассматриваемого метода. Способы выбора преобра­
зующей формы бyдyr рассмотрены ниже.
В
результате
интегрирования
системы
определяется
соответ-
ствующее значение 0/ Функции Q;
б) полученные значения Qi суммируются для определения оцен­
ки искомойвеiJиЧины n с весовыми коэффициентами pi:
(4)
Весовые коэффициентыI в соотношении
(4)
чтобы полученная оценка среднего значения
т.
е.
математическое
значением
ожидание
оценки
необходимы для того,
n
была несмещенной,
совпадало
бы
со
средним
П. Ниже по:казано, что соответствуюIЦИЙ выбор весовых
коэФФициентов в
(4)
позволяет получить несмещенную оценку при
любом чис.п:е испытаний и любой Функции О.
Пошем,
что
для
этого
весовые
коэффициентЫ·
должны
опредешiтьcя ·соотношением
где
нормировочный коэффициент, определяемый из условия ра­
Q -
венства единице среднего
штрих
(')
значения весовых коэффициентов
Q будет рассмотрен ниже. РаВенство р = 1
щmяется необходимым условием несмещенности оценки
вие легко получить, представив, что
Q. == Qi
для доказательства соотношения
в
=1,
обозначает транспонирование. Вопрос о вычислении норми­
ровочного коэффициента
случай,
р
котором
белый
шум
(5)
(4).
Эro усло­
== 1.
рассмотрим конечно-мерный
представлен
в
виде
ступенчатой
Функции
lЗЗ
=1, ... ,М,
j
где
't = Т/М -
чины ~ j
-
ширина одной ступеньки, М
число ступенек, а вели­
-
независимые случайные векторы, состоящие :каждый из
k
независимых случайных величин, распределенных по нормальному за-
кону с нулевым средним и дисперсией
1/'t.
распределен
с
по
нормальному
закону
-=
Вектор ~
IШотностью
{~i,~2, ... ,~M}
,
распределения
p;(~):
(6)
где
R; -
корреляционная матрица
вектора
~
размера
kM
х
kM,
равная
J
о
о
't
=
R;
1
О
О
't
О
1
О
't
1
-+
единичная матрица размера
kx k.
Компоненты вектора ~ = {~i,~2, ... ,~M}', являющегося ступен­
чатым приближением
11M(t) функции 11(t), определяются соотноше­
нием
м
~j
=
~j + 11j'
11j
= LKjm~m,
m=l
где К jm - элемент матрицы К = {К jm} , определяемый по формуле
K jm =
~
f f K(a,~)df3 da.
jt
[
m'!:
]
(j-l)t (m-l)'t
Orcюда
следует,
что
вектор
~
распределен
закону с IШотностью распределения Рс. (~):
134
по
нормальному
(7)
Rr.
где :корреляционная матрица
:вeкropa ~ определяется coorношением
~ =R; +.!'t (КК' + К' + К).
Условие
несмещенности 'оценки . Q
следующим образом: с одной стороны,
поСле
Q
замены
~
на
и
l;
.Q
I
может
бытъ
получено
= ~)p; (~M, с дрyroй
введения
весовых
-
коэффициентов
=I n(~)PI;(~~(~(~)~ . Заменяя в первом из этих интегралов
обозначение
получившееся
переменной
интегрирования
выражение
со
вторым
с
~
на
интегралом,
l;
и
сравнивая
получаем
условие
несмещенности оценки .Q в виде
Подставляя в это равенство плотности вероятности
(6)
и
(7),
получаем, что весовой коэффициент долж~н равняться
~~) =Q ехр(-i[~#{Я;-1 - Rr.-l}~] =
=Qexp(-~['t~'{К' + К + K~}~п =
=Qex~-~111'11 + 211,~п,
rде Q
1
=
1
~I2~Г2
-" веричина, не зависящая от ~ .
Возвращastсь к функu.иям:
'I1M(t)
и ~M(t), получаем
Устремляя Мк бесконечности, получаем соотношение
(5).
В заключение отметим, что в тех случаях, когда определение НОР­
мировочноro множителя
ero
-
1",,·
оценку Q ~ N
Q
затрудниreльно, в
(5)
можно использовать
N
L
;=1
р'
. При этом выражение (4) :Принимает вид
'
135
N
LOipi
n=
=....;~"":-_ _
.;:..1
Lpl
j=1
Способ...
3.
в...бора
преобразующей
форм....
Известно,
что
дисперсия оценки при использовании весового метода ~онте-Карло
сильно зависит от выбора преобразования для плотности веРОЯТНОС'i'и
[4].
В рассматриваемом методе она определяется способом выбора
преобразующей формы
(3).
В общем случае преобразующую форму
следует выбирать из условия минимума дисперсии получаемой оценки
Do
.. - _)2
о
~ inf,
=(о' р'
(8)
поскольку именно величина определяет, сколько реализаций потре­
буется для вычисления оценки с заданной точностью. В самом деле,
относительная погрешность оценки, полученной в результате выпол-
D
нения N
уменьшая
реализаций, равна М1::;: _ ~, и уменьшить ее можно,
n"N
DQ
и увеличиваяN. Последнее, однако, не вcerдa возмож-
но. Например, в задачах оценки малых вероятноетей,потребное число
реализаций
N
может оказаться слишком большим, и следует попы­
таться уменьшить величину п о . Следует отметить, что минимум в
(8)
не Bcerдa может быть найден. Поэтому наряду с точными решениями
(8)
следует также рассматривать эмnирические- алгоритмы задания пре­
образующей формы. Ниже npиведены
примеры, в одном из которых
используется оптимальное преобразование возмущений, в другом же
используется эмпирический подход, и за счет удачного выбора преоб­
разующей формы оказывается возможным существенно снизить по­
требное число реализаций.
4. Случай линейной систем... с посТоянн",ми коэффициентами и
фиксированн",м
временем.
Приводимый
ниже
пример
иллюстрирует
возможность выбора преобразующей формы в соответствии с усло­
вием
(8).
Рассматривается динамическая система, описываемая линейным
дифференциальным уравнением:
x(t) = A(t)x(t) + g(t~(t),
}
7I
k
х(О) = О, x(t) Е R , ~(t)E R ,
t Е [О,Т],
rдe
136
A(t)
и
g(t) -
т =
const,
матрицы соответствующей размерности.
(9)
n(Х(Т») = 1, Ф{х(Т») > Ф.,
.
Ф(х)
О, Ф{~(Т») s ф.,
в интегральной форме уравнение
(9)
= h'x.
имеет вид
f
X(t) = JW(t,t)g(t~('t}dt,
(10)
О
где
W(t,'t) -
решение матричного дифференциального уравнения
aW(t,r) ( ) ( )
at
W('t, 't) = 1.
-~~=АtWt,'t,
Отсюда следует, что
т
f
Ф = q'(t)~(t)dt,
(11)
О
где
q(t) -
импульсная переходная функция:
q(t) = h'W(T,t)g(t).
Из соотношения
(11)
следует, что
Ф
(12)
распределено по нормаль­
ному закону с нулевым средним и дисперсией 1~112:
Здесь и далее под нормой I~II функции q понимается норма
1
IИ ~ [I q'(t)q(t)dt):2.
Остановимся на вопросе выбора преобразующей формы. При лю­
бой реализации ~(t) добавок
rt(t)
в
(3) может, быть представлен в виде
1 т
J
p(t) = rt(t) - a,q(t). Из
.
'
м о
следует, что Слага.ем~ p(t) не будет влиять на распределение Ф при
rt(t) = a,q(t) + p(t) , где а, = - 2 q'(t)rt(t)dt ,
(11)
использовании весового метода, однако будет увеличивать дисперсию
весовых коэффициентов, откуда следует, что для оптимального преоб­
разования возмущений в соответствии с условием
(8)
p(t)
должно рав-
137
нятъся нулю. HeblИoro более сложными рассуждениями, которые здесь
не приводятся, можно ПОI<ВЗaТЬ, чro матричная ФУНICЦИЯ
влетворяющая условию
(8),
с точностью до :множителя
cr
K(t,e),
удо-
дomкнa иметь
цц
K(t,e) =cr q(t)q'(e) .
~~12
Исследуем зависимость дисперсии оценки от параметра
МатематичеСкое
ожидание
вероятность события Ф
> Ф. ,
-0 = -1-
П,
представляющее
очевlЩНО, не зависит от
ао
,,2
е
cr
собой
и равно
R2
I --dxlfl$---e
1 1 --
,fj;; R
где R =
оценки
cr.
2
,fj;; R
2
'
;i .
После замены ~(t) на l;(t) в соответствии с (3) имеем
т
I
т
о
о
J
ф,," q'(t)l;(t)dt "" (1 + cr) q'(t)s(t)dt ,
причем плотность распределенИя величины Ф теперь определяется со­
отношением
Дисперсия оценки
(13)
График зависимости функции
,
ш 2 (R cr)
S(cr,R) = Q(R)
от cr при различ-
ных значениях R приведен на рис. 1. эта функция показывает, на­
сколько может быть уменьшено число реализаций, потребное для
оценки П, так как это число пропорционально
N S(cr,R).
Ifl$
О
138
S(cr,R):
s (б, R)
10+
10 :
1
J
5
Ij.
R=J
10-2
5
If
'1.5
;i/i.
5
5,5
10-8
1.
Рис.
З3JIисимоcrь функции
S(O', R)
R
от о' при различных
значениях
Для обычноro метода Монте-Карло
функция
S(cr,R)
1
N == =.
= О,
При
n
равна единице и убывает с ростом а. При
неограниченно
минимально,
cr
возрастает.
можно
Выбирая
значительно
cr
так,
сократить
чтобы
cr ~
cr
00
S(cr,R)
потребное
=О
она
было
количество
реализаций.
Так, для
n
=:
2,8 .10-7
(R = 5)
возможно, выбирая cr
=:
4, умень-
шить потребное число реализаций примерно в 105 раз. Как видно из
1, функция. S(cr,R) в окрестности минимума меняется слабо,
рис.
поэтому точное
Отметим,
ero
определение не обязательно.
что
полученные
в
данном
справедливы для любой линейной системы вида
примере
(9)
результаты
и что для такой
системы оценка вероятности может быть вычислена roраздо
более
простым образом, например через корреляционную матрицу. Однако
предлатаемый метод направлен в первую очередь на решение нелиней­
HbIX задач, для которых данный линейный пример позволяет оценить,
какой экономии в' количестве реализаций можно достичь при опти­
мальном выборе преобразующей формы.
нелинейных задач
или
при
Следует ожидать,
неоптимальном выборе
что для
преобразующей
формы экономия в количестве реализаций будет не столь заметна.
5.
Пример 'примевевИJI формирующеro фильтра ДЛЯ заданИЯ преоб­
разующей
оценки
му
(8)
формы.
Точное
.решение
задачи минимизации дисперсии
во многих случаях является' очень сложной задачей, поэто­
наряду с точным ее решением следует рассмотреть также возмож­
ность выбора преобразующей формы эмпирическим ~:м. Пример та­
KOro подхода представлен ниже.
139
Дисперсия оценки
зависит, с одной стороны, от разброса
(8)
значений оцениваемой :величи:ны:, с дрyroй
- от дисперсии весовых
коэффициентов. Поскольку доля реализаций, для которых л,i '1: О, .
невелика, уменьшить дисперсию оценки (8) можно, y:вeJIИЧИВ3Я долю
таких реализаций. Orcюда следует, что npeобразующую форму нужно
выбрать
так,
системы
(1).
чтобы
увеличить
отклонения
фазовых
переменных
Эroго можно ДОСТИЧЬ, изменяя спеtcrpальную плотность
исходных возмущений.
Известно, что на отклонения фазовых переменных динамической
системы ТШIа
исходных
(1)
по-разному влияют участки спектральной плотности
возмущений.
Так,
ее
увеличение
в
области
достаточно
высоких частот пракгически не будет влиять на отклонения фазовых
переменных,
но
в
то
же
время
приведет
к
увеличению
разброса
весовых коэффициентов, что в сумме приведет к увеличению диспер­
сии оценки
(8),
что нежелательно. Поэтому при выборе npeобразу­
ющей формы следует увеличивать спектральную плотность возмущений
в
тех
областях,
которые
наиболее
сильно
влияют
на
отклонения
фазовых переменных динамической системы. Подобные раССy.JЩения
не
являются
полезными,
строгими,
так
:как
однако
MOIYТ
нелинейным системам
в
некоторых
применяться
случаях
:как
к
MOгyr
оказаться
линейным,
так
и
(1) .
. Целенаправленная
деформация спектральной плотности может
осуществляться, например, формирующим фильтром
it i(t) = F(~)"i (t) + G(t)~i (t),}
,,' (О) = "0'
(14)
что соответствует частному случаю преобразования вида
(3)
t
J
"/(t) = K(t, е) ~I (e~,
о
где
K(t,e)
удовлетворяет ураВнениям
aK(t е)
д/
F(t)K(t,e),
к(е,е) = а(е).
=
Orметим также, что в данном случае преобразованиевозмущений
осуществляется
О~е~t ~Т
,
не
на
квадрате
[О, Т] х [О, Т] , а
на
треугольнике
что более удобно для практической реализации статисти­
ческого моделирования, так :как для пр~образования возмущений не
требуется знать «будущих. (е
В
качестве
примера
> t)
возмущений.
применения
такого
подхода
рассмотрим
самый простой случай уравнения первого порядка. Требуется оnpеде-
140
лить верояmость превышения величиной х(Т) уровня Х•• Рассмотрим
дифференциальное уравнение
x(t) = -x(t) + .fi~(t),}
(15)
t е [-ао,Т1
описывающее стационарный случайный процесс. Orметим, что линей­
ность уравнения
(15)
в данном примере используется только для на­
хождения точноro решения и сравнения еro с результатами, получен-.
ными с помощью предлагаемоro метода:
в соответствии с приведенными выше рассуждениями
для
,,(t)
был ВЗЯТ формирующий фильтр вида
,,(t) = -,,(t) + d;(t),
t;::: О,
,,(0) = о.
с помощью параметра
d
можно бьmо определять степень дефор­
мации спектральной плотности. Увеличение этоro цараметра равно­
мерно увеличивает спектральную плотносТь в . области низких частот,
не изменяя ее в области высоких частот . Усиливать высокочастотные
составляющие возмущений, очевидно, не имеет смысла, так как систе­
ма
(15)
к ним не чувствительна.
40 000
Были выполнены серии испытаний по
личных значений параметра
d
при Т
=
10
реализаций для раз­
(выбор именно такоro зна­
чения Т объясняется тем, что для практических задач динамики поле­
та характерное отношение времени Т в
равно
z 10 + 30).
(1)
к времени
автокорреляции
Полученные результаты представлены на рис.
2-5 .
••••••••••• ••
06 а
"i'I" Вее
• • • ••••
•••• •
511
"QD
10* .о
о
5
В
В D
О е е о
•
10'
D
D
о
о
Рис.
о
• -
учета
2
веСОВОЙ метод Моиre-Карло,
40 000
реализаций,
d" 2,0;
вероятнocrь реализаций в процессе моделирования без
весовых
коэффициеJПОВ;
О
-
Иcrинноо
значение
вероятности
141
•••••••• ••• •• • • • •
08.
888
8.
8
о
1Ii'
• •••••
911а
986
Q
•
8
D
о
а
о
D
Рис.
3
11 е
Обозначении тахие же, хак на рис.
2
II(СНliI6IИ6. 4tfН/RmHOem6
200 1
10-2
"г·
10-/
10-111
10-1
10-11.
о
о
о
о
о
Рис.
4.
Оценка величины ornосительной ПОJpeшности
ДЛJI различных значений парам:етра
На рис. 2 и
при значенияхd
3
d
представлены оценки вероятностей, полученные
и 2,0 соответственно. По оси Х отложено ис-
=1,5
тинное значение вероятности, вычисляемое :как
По оси У отложена оценка вероятности, полученная с помощью
предлагаемого метода при соответствующем значении naраметра
d,
а
также истинное значение вероятностИ и вероятность реализации собы­
тия х(Т) ~ х. в процессе моделирования. для значения
совпадение оценки с
d =1,5
истинным. значением наблюдается
хорошее
вплоть до
значений Р ~ 10-9, а для d =2,0 - вплоть до значений Р ~ 10-11. На
рис. 4 показана оценка величины относительной по:rpeшности, прихо­
дящейся на одно испытание, для различных значений парам:етра d:
·е·
-1
.р
142
-.
;
pi,
х ~ х.,
О,
xi<x..
lO-,z
1051
0/(tно60еион 'еронmносmll
10lf
10-8
10-1
10-711
А
tl.=l,J
1,5
...
1,7
Q
2
J
D
о
Рис.
5.
1012
Оценка потребноro числа реалиэаций для раз­
личных значений параметра d
Учитывая, что с ростом числа реализаций поrpeшность оценки
вероятности убывает как cr ~ а/Ш, потребное число реализаций может быть оценено :как N ~ а 2 . Оценка этой величины для различных
значений параметра d приведена на рис. 5. Из него следует, что даже
для значений вероятности Р ~ 10-11 потребное число реализаций со­
ставляет N ~ 105, в то время как при применении обычного метода
Монте-Карло потребовалось бы порядка 1011 реализаций. Таким об­
разом, для данного примера предлагаемый метод позволяет сократить
число испытаний на пять порядков.
Работа выполнена при поддержке РФФИ.
ЛИТЕРАТУРА
1.
С т р а т о н о в и ч Р. л. Новая форма записи стохастических юrге­
гралов и уравнений
//
Вестник МТУ. Серия мат. и мех.-
1964, N.! 1.
Кузьмин В. п.,
Ярошевский В. А. Асимптотические
оценки вероятностей больших случайных отклонений фазовых координат ди­
2.
намических систем от средних значений // Ученые записки ЦАТИ.- 1990.
Т. XXI, N.! з.
з. Д ы н н и к о в А. А. Асимптотическая оценка вероятностей больших
случайных отклонений фазовых координат динамической системы при' нали­
чии фазовых ограничений // Ученые записки ЦАТИ.- 1994. Т. ХХУ, N.! 4.
4. М и хай л о в г. А. Оптимизация весовых методов Монте-Карло.­
М.: Наука.- 1981.
Рукопись поступила
27/IX 1994 г.
1/--страниц
Пожаловаться на содержимое документа