Использование метода существенной выборки при моделировании автоматической посадки Текст научной статьи по специальности «Математика»

Том XXX V

УЧЕНЫЕ ЗАПИСКИ ЦАГИ 20 0 4

№ 3 — 4

УДК 629.735.33.015.077

ИСПОЛЬЗОВАНИЕ МЕТОДА СУЩЕСТВЕННОЙ ВЫБОРКИ ПРИ МОДЕЛИРОВАНИИ АВТОМАТИЧЕСКОЙ ПОСАДКИ

А. А. КУЗНЕЦОВ, В. П. КУЗЬМИН

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

1. Оптимальное распределение для оценки вероятностей. Оценка вероятности методом Монте-Карло сводится к определению многомерных интегралов. Пусть задана некоторая функция у(х), где х = х1, хк — вектор случайных величин с плотностью распределения

вероятностей/(х), тогда вероятность Р(у > у ) сводится к вычислению интеграла

Р( У > У *) = | Ф( У) / (*)«Х = М [ф],

|0 1 бе у < у*,

где О — область изменения переменных х1, ..., хк, ф(у) = < М[ф] —математическое

[1 1 бе у > у*,

ожидание случайной величины ф.

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

Оценка вероятности по N реализациям вычисляется по формуле для среднего арифметического:

N

р=--------------, (1)

N

где X — независимые случайные реализации, моделируемые по заданной плотности распределения /(х).

Дисперсия оценки, очевидно, определяется как В[Р^ = В0, Где В0 = В[ф] = Р0 - Р02 —

дисперсия на одно испытание, а Р0 — истинное значение искомой вероятности. Поскольку функция ф принимает значения 0 и 1, то расчет вероятности методом Монте-Карло (1) сводится просто к вычислению относительного числа реализаций, в которых выполнилось условие у > у .

В качестве меры точности оценки вероятности методом Монте-Карло будем использовать отношение среднего квадратического отклонения оценки к ее среднему значению. Задав

требуемую относительную точность

В

Р

N >

< в, легко получить потребное число реализаций Вп

Ро2в2

Основная идея, на которую опираются методы уменьшения дисперсии оценки, заключается в замене исходной плотности распределения вероятностей /(х) на другую функцию £(х). Различные практические методы уменьшения дисперсии оценки методом Монте-Карло отличаются

в основном лишь способом подбора подходящей функции £(х) [1]. Вероятность Р(у > у) при такой замене может быть переписана в виде:

р (у >/)=/

ф( У)

/ (х) g (х)

g (х)ёх.

Оценка вероятности методом Монте-Карло определяется суммой

^ У (х х)

g (хх)

Р =

N

£ф[ У( х х)

х=1

N

где X — независимые реализации случайного вектора х с плотностью распределения £(х). Дисперсия оценки на одно испытание составит

В1 = |ф[ У(х)]

(2)

Основная цель метода существенной выборки — выбрать новое распределение вероятностей так, чтобы сделать дисперсию оценки минимальной.

Пусть все возмущения разделяются на две группы х = (х1, х2). Эти группы отличаются тем, что мы хотим изменять распределение вероятности У1(х1) только для вектора х1. Плотность совместного распределения вероятностей представляется в виде /(х1, х2) = У1(х1) /(х2|х1).

Найдем новое распределение g(x1), обеспечивающее минимум дисперсии оценки. При замене распределения /1(х1) на g(x1) оценка вероятности находится по формуле

N

/1 ( х/ )

Р = ^ф[У(x1,х2)] / Л • х=1 g ( х1 )

Дисперсия на одно испытание составит

В1 =Цф2(У)/2 (х2 1 х1)

/12 (х) g(х1)

dx1dx2.

(3)

Учитывая, что ф2(у) = ф(у), а |ф(у)/2 (х2 | X!)ёх2 = Р(у > у*|х^ — условная вероятность

*

неравенства у > у при заданном значении хь получим

D = {Р(у > /к)^-TYdxi -P2. J g (xi)

Функция g(x1), которая доставляет минимум дисперсии (3) при условии {g(xi)dxi = 1, определяется соотношением

g (Xi)=_ЯЩйлх^_. (4)

JVP(У > у Vi) f (xi)dxi

Дисперсия оценки на одно испытание (3) при этом примет следующий вид:

D = {V Р(У > y*|xi) f (Xi)dxi - P02. (5)

Рассмотрим два простых случая.

1. Функция y(xi, х2) зависит только от вектора х2, т. е. не зависит от переменных, для которых изменится распределение вероятности. Тогда Р(у > у |xi) = const = P0 и дисперсия на

-

одно испытание составляет Di = Р() - Р() = D00, т. е. равняется исходной дисперсии.

2. Функция у(хь х2) зависит только от вектора xi. Это предположение может соответствовать тому, что мы меняем распределение вероятностей для всех переменных. Тогда

ш . f0, ааа у(xi) < у*,

Р( У > /|xi) = 1 ф

[i, ааа у(xi) > у ,

{VР(У > /Vi)f (xi)dxi

= Ро.

2

и дисперсия на одно испытание принимает наименьшее значение = Р0 - Р0 = 0.

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

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

2. Оценки больших отклонений вертикальной скорости и дальности при автоматической посадке. Описанная выше методика применяется для оценки распределений вертикальной скорости и дальности в момент касания ВПП при автоматической посадке самолета.

Рассматривается упрощенная модель автоматической посадки самолета в вертикальной плоскости при действии ветровых возмущений.

Различные математические модели ветра в приземном слое атмосферы (Н < 300 м) отличаются параметрами и видом спектральной плотности для порывов ветра, но имеют некоторые общие свойства [2]. Ветер представляется в виде суммы двух слагаемых: зависящего от высоты горизонтального градиентного ветра и(Н) и произвольно направленного турбулентного ветра

с быстрой сменой значения и направления. Градиентный ветер задается независимыми проекциями на оси земной системы координат:

их (н) = пхф(И), иг (И) = иг ф( И).

Общепринятым является логарифмический профиль ветра в зависимости от высоты в приземном слое

где Н — высота центра масс самолета в метрах, Н0 — параметр шероховатости подстилающей поверхности с характерными значениями от 0.05 до 0.1 м в зависимости от рельефа местности вблизи аэродрома.

Случайные величины их и иг являются независимыми с распределениями, которые могут быть аппроксимированы усеченными гауссовскими законами с параметрами [2]:

Три составляющие турбулентного ветра (порывы) считаются гауссовскими процессами со спектральной плотностью в форме Драйдена [2] и случайными для каждой реализации среднеквадратическими отклонениями, которые пропорциональны модулю градиентного ветра на высоте 10 м:

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

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

где ^1, ^2 и являются независимыми, нормально распределенными величинами с параметрами

(0.1). В формуле (7) первое слагаемое соответствует воздействию турбулентности, а второе — определяет влияние продольного ветра и, следовательно, переменная ^1 соответствует продольной составляющей градиентного ветра, ^2 — боковой составляющей, а ^3 — турбулентности. Турбулентная компонента задается с помощью одной переменной, несмотря на то, что турбулентный ветер задается большим числом случайных величин. Замена случайных величин одной объясняется тем, что эта величина является их линейной комбинацией, а случайные величины имеют нормальное распределение.

Для каждой из оцениваемых величин, дальности (Ь) и вертикальной скорости (V), строится своя линейная модель. Единственным отличием и единственным существенным параметром линейных моделей является коэффициент, стоящий перед вторым слагаемым в следующих формулах:

= 18 ( Н/Н 0 )

их =-2.7 1 и2 = 0 СТи =СТих =СТиг = 3.75 м/с.

(6)

(7)

(^1 + их/°и)2 +^2 + ЬЬ^1 = кЬ [^зи + ЬЬ^1 ].

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

проведении расчетов коэффициенты Ъу , Ъь задавались в виде Ъьу = аь у

Л

г- ^

Ux

\аы У

. В этом

случае коэффициенты ау , аь с точностью до знака совпадают с отношением средних

квадратических отклонений соответствующей координаты, вызванных только градиентным ветром и только турбулентностью.

Рассмотрим случай, когда изменяется распределение вероятностей только модуля ветра и. Так как параметры их и иг являются независимыми, то функция совместного распределения компонент градиентного ветра представляется в виде произведения двух функций распределения

/(их) и /(иг):

f (Ux, Uz ) = f (Ux )f (Uz ) =-

1

2ла„

exp

2 а2

іux - ux )2

Для моделирования существенной выборки по одной переменной, модулю ветра, необходимо перейти к новым переменным

Гux = u cos ф,

\uz = u sin ф.

Тогда функция совместного распределения представится в виде:

f (u, ф) =

1

2па:

exp

1

2а2

" 2 2 iu cos ф-йх) +і u sin ф)

>u.

Плотность распределения модуля ветра /(и) есть интеграл от совместной плотности распределения /(и, ф):

—2 2 1 ux: 1 u

2п 1 --^ -Л^_ 2п

f (u) = f f (u, фМф = —- e 2 а“ e 2 а“ u Г euux cos фdф, J 2а2 J

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

f (ф I u) =

f (u, ф) f (u)

pUUxcosф

2n

f euUxcos фd ф 0

Будем оценивать вероятности P(R > R*) для безразмерных переменных Rl =

L - L

Of

и

V - V

R = У У RVy =-----

а

. Для безразмерных переменных модель (7) имеет вид

УУ

R = ^з

г + a {u cos ф- ux ).

д/2 + ^xlau )2

Условное распределение вероятностей P(R > R* u) определяется очевидным образом

u

Рис. 1. Оптимальные плотности распределения модуля Рис. 2. Оценки эффективности существенной выборки по ветра модулю ветра

2п

Р(Я > Я*\и) = | Ф

Я* _а(и008ф-и

/ (ф\и) Л ф,

(9)

ТО

1 С _ 2 2

где Ф( х) = .— I е 'Ж — интеграл Лапласа.

-\/2п

X

Новая плотность распределения вероятностей g(u) определяется соотношением, аналогичным (4):

^и) = Р(Я > Я* \ и)/(и).

(10)

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

N0 = ЩР]

А[Р ].

N

(11)

На рис. 1 приведены характерные распределения вероятностей модуля ветра g(u) в зависимости от безразмерного уровня отклонения фазовой координаты. Для сравнения, на этом же рисунке приведена исходная плотность распределения /(и).

Рис. 2 иллюстрирует изменение потребного числа реализаций в зависимости от величины оцениваемой вероятности и параметра а.

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

В табл. 1 приведены значения относительных погрешностей оценок вероятностей для различных уровней при условии, что все оценки получаются с помощью одного оптимального распределения вероятностей, вычисленного для уровня Я = 6.

Т аблица 1

Относительная точность оценок для разных уровней е(Я )/б(6)

Я* 1 2 3 4 5 6 7 8 9

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

а = _0.5, Р0 = 1.5-10-4 14.2 10.3 3.6 0.96 0.72 1.0 1.7 2.9 5.5

а = 0, Р0 = 9.2-10-5 7.2 6.2 2.6 0.84 0.69 1.0 1.7 3.0 5.7

а = 0.5, Р0 = 1.6 10-5 3.5 4.6 2.7 0.82 0.60 1.0 2.0 4.1 9.0

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

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

Рассмотрим изменения в определении новой плотности вероятностей, которые позволят устранить данный недостаток. Пусть требуется построить распределение вероятностей для диапазона уровней [0 — Я ] и задана относительная точность оценок в зависимости от уровня в(Я*). Распределение потребного числа реализаций для заданного уровня будет иметь вид

gN (и, Я*) = g(и,Я*) 2 ^Р2 т . (12)

Р02(Я )в2(Я )

Такие плотности приведены на рис. 3 для точности оценок, соответствующих значению 1 = 5.

82( Я)

Для того чтобы обеспечить требуемую точность одновременно для всех уровней из диапазона [0 — Ятах], достаточно задать плотность числа реализаций, соответствующую огибающей семейства плотностей (12):

gl(u,Ятах) = . тах gN (u,Я*).

Я =[0,Ятах]

Потребное число реализаций, очевидно, составит N = | g1(u, Ятах)Ли.

В табл. 2 приведена величина (N1 _ N)/N, характеризующая увеличение потребного числа реализаций при использовании огибающей семейства плотностей для диапазона [0 — Ятах] по сравнению с использованием плотности для одного уровня Ятах.

Таблица 2

Относительное увеличение потребного числа реализаций при использовании огибающего распределения (%)

-30 -20 -10 0 10 20

Рис. 5. Оптимальные плотности распределения компонент ветра их и щ

^шах 2 3 4 5 6 7

а = -0.5 5.3 4.2 3.2 2.3 1.6 1.1

а = 0 4.0 3.1 2.2 1.6 1.1 0.8

а = 0.5 2.3 1.6 1.0 0.7 0.4 0.2

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

На рис. 4 иллюстрируется эффективность метода для оценки вероятности при различных значениях параметра а. Здесь приведено отношение потребного числа реализаций для метода Монте-Карло с исходными и оптимальными распределениями вероятностей при одинаковой точности оценок для различных значений оцениваемой вероятности. При значении а = 0 эффективность выборки по двум переменным совпадает с эффективностью выборки только по модулю ветра, поскольку в этом случае модель (8) зависит только от модуля ветра. При значении а Ф 0 эффективность метода по двум переменным несколько выше, в первую очередь — при а > 0

На рис. 5 приведен пример оптимального распределения вероятностей. Основной особенностью оптимальных распределений является то, что случайные переменные их и иг становятся зависимыми, тогда как для исходного распределения эти переменные являются независимыми.

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

Т аблица 3

Средние значения и средние квадратические отклонения

Ы[Ь], м М[Уу], м/с аь, м Сту,, м/с

Влияние турбулентности 418.3 -0.73 18.9 0.14

Влияние всех возмущений 422.8 -0.84 45.58 0.33

Т аблица 4

Зависимость параметров приземления от величины продольного ветра

их, м/с -8.1 -5.4 -2.7 0 2.7 5.4

Уу, м/с -0.84 -0.78 -0.70 -0.64 -0.61 -0.57

На основании этих результатов были определены приближенные значения параметров аь « 0.6 и ау « -0.4. Для получения оценок параметров распределений были выполнены расчеты

по стандартному методу Монте-Карло объемом N = 105 реализаций и с использованием новых распределений по двум переменным. При проведении статистических расчетов задавались уровни с интервалом АЯ = 0.5 и использовались оценки для тех уровней, которые были превышены не менее, чем в пяти реализациях.

При использовании метода существенной выборки определялось распределение вероятности, обеспечивающее получение оценок для всего диапазона уровней [0 — Яшах], а величина уровня Яшах подбиралась так, что потребное число реализаций составляло N1 « 105. Таким образом, для оценки эффективности метода существенной выборки сравнивались минимальные значения найденных вероятностей.

На рис. 6, 7 приведено сравнение распределений вероятности для вертикальной скорости и дальности, полученных стандартным методом Монте-Карло (сплошные линии) и с использованием измененных распределений вероятностей (кружки). Данные результаты показывают хорошее совпадение оценок, полученных с использованием метода существенной выборки и стандартного метода Монте-Карло.

Сравнение предельных вероятностей, полученных двумя методами, показывает, что использование метода существенной выборки позволяет сократить объем расчетов более, чем в 100 раз.

Рис. 6. Вероятности превышения характерных значений дальности

Рис. 7. Вероятности превышения характерных значений вертикальной скорости

Оценки, полученные ранее для модельной задачи (см. рис. 4), для значений вероятностей (~10-6 г 10-7) дают сокращение объема расчетов также примерно в 100 раз. Совпадение оценки эффективности метода для модельной задачи с полученной в результате статистического моделирования дает основание полагать, что вычисленные с помощью модельной задачи распределения вероятностей составляющих градиентного ветра оказываются близкими к оптимальным.

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

ЛИТЕРАТУРА

1. Ермаков С. М. Метод Монте-Карло и смежные вопросы.— М: Наука.— 1975.

2. Кузьмин В. П., Ярошевский В. А. Оценка предельных отклонений фазовых координат динамической системы при случайных возмущениях. — М.: Наука.— 1995.

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