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

Том XXIX

19 98

№1-2

УДК 629.735.33.015.077

ОЦЕНКА ВЕРОЯТНОСТИ УСПЕШНОЙ ПОСАДКИ НЕМАНЕВРЕННОГО САМОЛЕТА ВЕСОВЫМ МЕТОДОМ МОНТЕ-КАРЛО

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

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

В общем случае процесс посадки самолета может быть описан системой стохастических дифференциальных уравнений:

где х(г) — вектор фазовых координат динамической системы, включающий в себя как фазовые координаты собственно самолета и его си-

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

*(*) = /(*(0« 0 + £(*(0. х(г)єД", Ф)єЛк, / є [0, Г], • х(0) = х0,

(1)

стемы управления, так и формирующие фильтры, используемые для моделирования атмосферной турбулентности. Я", Лк — евклидовы пространства соответствующих размерностей. Через с(0 обозначен вектор некоррелированных гауссовских белых шумов с единичной интенсивностью, подаваемых на вход формирующих фильтров, g — матрица соответствующей размерности. Стохастическая система (1) понимается в смысле Стратоновича [1]. Момент Т окончания процесса посадки определяется из условия равенства нулю высоты полета. Формально это условие может быть представлено в виде ф(л,(7’)) = 0. Область допустимых значений для фазовых переменных в общем случае может быть задана неравенством Ф(х) < Ф„, откуда следует, что Рл п = Р{ф{х) > Ф*).

Отмстим, что эта вероятность тождественно равна среднему значению функции

|1. ФМ,Ф..

[О, Ф(х)<Ф*,

Р.1.П

Согласно существующим нормативам эта вероятность не должна превышать 10_6 +10-8. Оценка столь малых вероятностей требует вычисления большого числа реализаций, поскольку искомое событие в процессе моделирования встречается крайне редко. Для уменьшения требуемого числа реализаций предлагается использовать модификацию весового метода Монте-Карло, предложенную в работе [2]. Идея этого подхода состоит в том, чтобы, изменяя структуру возмущений и увеличивая частоту реализации искомого события в процессе моделирования, добиться уменьшения необходимого числа реализаций. При таком подходе искомая вероятность уже не будет равна частоте событий в процессе моделирования и для ее оценки необходимо использовать пересчет с помощью весовых коэффициентов. Этот подход был практически реализован для данного класса задач в [2] и может быть описан следующим образом.

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

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

(здесь и далее верхний индекс означает номер реализации), определяемый соотношением

Т„,

?'(/) = л'»+ £''(/), т/(0= |*М)£'(еуе, (2)

о

где А'(лб) — кусочно-непрерывная матричная функция, заданная на

квадрате [0, Тт] х [0, Тт], Тт — максимальное время интегрирования

системы (1). Функцию К(/,б) будем называть преобразующей формой.

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

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

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

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

шеупомянутое условие несмещенности оценки М[р\ = 1.

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

поскольку именно величина 2)п определяет количество реализаций, требуемое для вычисления оценки с заданной точностью. В самом деле,

(3)

/

р' = Qp' , р' = ехр

V

о

N

равно единице. Выражение (3) принимает

этого

выражения следует, что среднее значение оценки значению весовых коэффициентов

Л/Jq] равно среднему Отсюда вытекает вы-

2

(5)

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

шая Х>а и увеличивая N. Последнее, однако, не всегда возможно. Например, в рассматриваемой задаче оценки вероятности летного происшествия потребное число реализаций ТУ оказывается слишком большим, и следует попытаться уменьшить величину /)а. Следует отметить, что точный минимум в (5) может быть найден лишь в достаточно простых случаях, например для линейных задач. В то же время нахождение этого минимума в общем случае представляет собой задачу намного более сложную, чем исходная задача оценки вероятности, поэтому следует также рассматривать эмпирические алгоритмы задания преобразующей формы. Именно этот подход использован ниже для оценки вероятности Рл п. Как было отмечено выше, уменьшить дисперсию оценки (5) можно за счет увеличения доли реализаций, в которых имеет место интересующее нас событие. Этого можно достичь, например, изменяя спектральную плотность исходных возмущений так, чтобы увеличить дисперсию фазовых переменных в момент окончания процесса посадки.

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

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

Для уменьшения дисперсии оценок искомых вероятностей использовалось преобразование возмущений с помощью формирующего фильтра первого порядка, т. е. в уравнение (1) вместо вектора белого шума £(/) подставлялся вектор ¿¡(/), определяемый как С(0 = 5(*) + л(0» где г|(/) — решение дифференциального уравнения

Выбор такой формы преобразования объясняется тем, что в данной задаче на отклонения фазовых переменных влияют лишь низко-

и уменьшить ее можно, умень-

(6)

частотные составляющие белого шума. Фильтр (6) увеличивает низкочастотные составляющие возмущений. Степень увеличения определяется параметрами ¿их. Ниже приведены результаты моделирования, полученные при различных значениях параметра <1. Постоянная времени, равная т = 12 с, была выбрана экспериментально на основе результатов, полученных при моделировании для различных значений параметра т, в ходе которых оценивалась степень влияния различных частот спектра возмущений на отклонения фазовых переменных, и в дальнейшем не варьировалась. Выбранная величина была сопоставима с характерной длительностью моделируемого процесса посадки (~20 с).

Для различных значений параметра сI (¿/=0,6; 1,2; 1,8; 2,4) были выполнены серии по 105 испытаний, в которых определялись оценки вероятностей таких событий, как посадка со слишком малой (1 < ¿т|п, Хтт — минимально допустимое для посадки расстояние от среза ВПП) и слишком большой (Ь > £тах, £тах — максимально допустимое для посадки расстояние от среза ВПП) дальностью, посадка с нарушением ограничения на вертикальную скорость и т. д., а также дисперсии этих оценок. В качестве примера на рис. 1—3 представлены результаты, полученные для случая посадки со слишком малой дальностью. На рис. 1 представлены оценки вероятности посадки с нарушением ограничения

Рис. 1

по минимальной дальности от среза ВПП (номинальное значение для этой величины составляет около 400 м), выполненные при различных значениях параметра с!. По оси абсцисс отложено пороговое значение дальности Ьтт, по оси ординат — вероятность посадки с дальностью меньше, чем Хт;п.

Отметим, что на рис. 1 и 2 представлены все полученные оценки вероятности, в том числе и те, для которых в процессе моделирования встретилась всего одна реализация. Очевидно, что такие оценки не являются достоверными: для получения достоверных оценок требуется десять и более реализаций. Тем не менее видно, что применение весового метода позволяет получить оценки в более широком диапазоне вероятностей, чем обычный метод Монте-Карло. Как и следовало ожидать, полученные при различных значениях параметра с/ оценки совпадают (в области, где эти оценки достоверны), хотя частота посадок с недолетом при моделировании существенно увеличивается (рис. 2) и зависит от с/. Значения отношения дисперсии оценки к дисперсии, соответствующей'обычному методу Монте-Карло, представлены на рис. 3. По горизонтальной оси отложена оценка вероятности Рл п, а по вертикальной — отношение дисперсий:

^ = Ю5. (7)

л.п/-** /

Функция 5(РЛ П) позволяет оценить, насколько может быть уменьшено число реализаций, необходимых для оценки вероятности Рл п, по сравнению с обычным методом Монте-Карло. Так, для вероятности Рл п ~ 10-6 из рис. 3 значение ¿'(/>л.п) составляет 8 ~ ДО-2, т. е. потребное число реализаций может быть сокращено в 100 раз. Учитывая, что в зависимости от требуемой точности при применении стандартного метода Монте-Карло для оценки вероятности требуется

N * ^ ^, Р « 1 реализаций, получаем, что для оценки указанной

вероятности предлагаемым методом потребуется порядка N « 3 • 104 3 • 105 испытаний. При Рл п ~ 10-8 имеем 5 ~ 3 • 10'4 и соот-

ветственно потребное число реализаций N а Ю5 + 106, что существенно меньше, чем при применении стандартной процедуры Монте-Карло N * 3 • 108 3 • 109. В случае «большой» вероятности Р:[ и ~ 10~2 данный

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

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

Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (код проекта 95-01-01309).

ЛИТЕРАТУРА

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

2. Д ы н н и к о в А. А. Применение весового метода Монте-Карло к исследованию редких событий в динамических системах, возбуждаемых белым шумом // Ученые записки ЦАГИ. — 1996. Т. XXVII. № 3—4.

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

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

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