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

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

№ 1

УДК 629.735.33.015.073

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

В. П. Кузьмин

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

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

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

2. Основные соотношения для оценки вероятности. Рассмотрим стохастическую систему вида

§-/(*> + *. 0)

где х — я-мерный вектор фазовых координат, f(x) —векторная функция, | — r-мерный вектор некоррелированных гауссовских белых шумов единичной интенсивности, g—матрица соответствующей размерности.

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

Без потери общности можно считать, что начальным условием для уравнения (1) является х(0) =0.

Пусть условие нормального функционирования системы, описываемой уравнением (1), сводится к неравенству для одной, например, первой координаты

*1 (0-0-

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

Пусть х = 0 устойчивое положение равновесия, а время Т достаточно велико, т. е. Т^>Тзат, где Гзат — время затухания переходных процессов в системе (1) без случайных возмущений. Если для случайной функции Xi(t) В некоторый момент t* выполняется условие Xi(t*)=a, то данная реализация не рассматривается при t>t*. Для оценки вероятности Р(а) существенным является дифференцируемость в среднем квадратичном случайной функции Xj(/) [3]. Если функция Xi(t) является дифференцируемой, то оценка вероятности пересечения заданного уровня может быть получена при больших значениях а через среднее число выбросов в единицу времени N.

Среднее число выбросов для стационарного процесса определяется формулой [2]

00

N= I” л:«>2(x,x)dx, (2)

о

где (дг{х, х) —совместная плотность распределения вероятностей переменной х и ее производной. Вероятность пересечения заданного уровня за время Т может быть определена по формуле

P(a)=\-e~NT, которая справедлива при Л/Т<С 1. В этом случае

Р (а) ~ ~NT.

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

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

Рассмотрим теперь случай недифференцируемой составляющей х1у), описываемой одномерным уравнением типа (1)

Уравнение Фокера — Планка для плотности вероятности, а также граничные и начальные условия имеют вид:

(0,(0, ¿) = 0,

где 6 (х) — дельта-функция.

Если отсутствует ограничение на изменение переменной хх^), то уравнение (4) имеет стационарное решение

где А —нормировочный множитель.

Решение нестационарного уравнения (4) при больших значениях а и больших t может быть представлено асимптотическим выражением [2]

где Хо—параметр, Дш(х1) —отличие плотности вероятности от стационарной, обусловленное граничным условием (4). Подставляя (5) в уравнение (4) можно по аналогии с [2] получить

Отметим, что формула (6) получена в предположении существования стационарного решения уравнения (4) без ограничения х1(1)<а. При оценке вероятности пересечения заданного уровня поведение функции /1(х1) при Л'1 >а несущественно. Формула (6) будет справедлива при произвольном поведении функции ¡1(х) при х>а, при выполнении условия

Формулу (6) можно использовать и в случае многомерной системы (1), если ограничение задано для одной координаты. На этот факт можно смотреть как на результат приведения многомерной системы (1) к стохастически эквивалентной (3) [3].

Рассмотрим возможное упрощение формул (2) и (6). Предположим, что система (1) допускает линеаризацию при малых отклоне-

(3)

Ш ~ ~2

(4)

—СО

® (■*!» О К (■*])+ Д°> (*1)],

(5)

(6)

а

х

— 00

Из (5) можно получить

Р Яй 1 - е-ч

(7)

ниях х от положения равновесия. Тогда одномерная плотность распределения вероятностей может быть представлена в виде соо^) =

. к3 ио

= _е 2 > где ^(0) =0 и Я(Хх) линейна при малых Х\. В этом слу-

чае можно приближенно задать нормировочный множитель А и величину gi как для линейной системы:

где <зХ1 — среднее квадратическое отклонение, а т*, — интервал корреляции случайной функции х^) ¡[4].

С учетом сделанных предположений формула (2) для среднего числа пересечений заданного уровня дифференцируемой составляющей х1 (I) примет вид:

где (х)—одномерная плотность распределения вероятностей переменной хь а — среднее квадратическое отклонение производной Х\

Как видно из формул (2) и (6), для определения вероятности превышения случайным процессом xt (t) заданного уровня xt (t),= a необходимо определить одномерную плотность вероятности переменной Xj в окрестности значения х^а или, что тоже самое, следует определить зависимость R(xi) в окрестности значения х^ — а {см. (8), (9)].

При больших значениях а зависимость R(a) может быть определена для многомерной системы (1) из решения вариационной задачи [5]:

Отметим, что задача (10) может решаться на интервале времени [0, Т4], причем величину Т1 достаточно выбрать превышающей время затухания переходных процессов Т3ат в системе (1), при этом величина Т1 может быть значительно меньше, чем Т.

3. Численный метод определения плотности вероятности. Рассмотрим численный метод решения задачи (10), который позволяет также приближенно определить и величину интервала корреляции переменной %х и ее среднее квадратическое отклонение а*. Пусть для простоты 5(0 скалярная функция. Будем искать решение задачи (10) в виде

И*(а)

2

(8)

2 тс *

(9)

Т г

(10)

при условии

•* = /1 С*) + gt, *1 (0) = о, (Т) = а.

m

5 (о=2 (о

(и)

5—«Ученые 'записки» № 1

65

где <fk(t) — ортонормированйые функции, Причем <рi(0= ;

■¡л¡(¿) = xl(t) — реакция системы на единичный импульс Ц() = Ь(0).

Если система (1) линейна, то ni(f) —импульсная переходная функция, и решение задачи (10) имеет вид

£(0 = Доданное обстоятельство и определяет выбор функции <Pi(i)- Величина R определяется из решения алгебраического уравнения

хг(Т, R) = a.

Б общем случае задача (10) с учетом (11) преобразуется к задаче

ТП

нелинейного программирования /?=min^c|, при условии

°ь fe=l

^1 (¿1> • • ч Crm Т) = а.

Учитывая сделанное ранее предположение о возможности линеаризации системы в окрестности положения равновесия, получим приближенные формулы ДЛЯ определения величин Хх и Ох-

J (12)

г т

(М (¿і) Iа (¿1 + т) dty dz

О

4

J J

a

(13)

Точность решения вариационной задачи (10) будет зависить от удачного выбора системы функций ф1; <р2, ..срт-

Пусть задана некоторая функция |о(0- Соответствующее ей значение Хі(Т) удовлетворяет граничному условию

хх (Г) = а.

т

Обозначим | ЦШ = /?£.

о

Пусть |(0—случайная реализация белого шума. Отметим, что при численных расчетах §(0 задается в виде кусочно-постоянной функции

Е(,)=Йг

при ¿*<¿<¿»+1, і=1, А£— шаг разбиения отрезка [0, Г]; Сі — не-

зависимые случайные нормально распределенные величины с единичной дисперсией и нулевым средним.

где е — малый параметр. Если |0(/)—оптимальное решение задачи (10). то при малых значениях е в первом приближении должно выполняться условие

Если |0(0 не является экстремальной для задачи (10), то можно определить поправку к величине фужкционала [6]

где Ох = М {[л^ (^ (¿)) — а]2! —дисперсия величины хи обусловленная воздействием случайных возмущений 5;^) (14) при е = 1.

Величина йх*/дЯ0 определяется численно путем интегрирования системы (1) для реализации возмущений, заданных в виде £(0=#£о(0 для нескольких значений Я. Величина Их может быть определена методом статистических испытаний (Монте-Карло).

4. Пример. Рассмотрим применение описанной выше методики для определения вероятности превышения углом атаки самолета заданного значения.

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

где а — угол атаки, —угловая скорость тангажа, —составляющая ветра, перпендикулярная к скорости полета, 1 — вспомогательная

переменная, 8„ = ЛР 8В> бв — отклонение руля высоты, V — скорость полета, ¿ — масштаб турбулентности, аут — среднее квадратическое отклонение скорости ветра, Vа, /И®, Мш —постоянные коэффициенты, определяемые аэродинамическими характеристиками самолета,

7

о

т

о

ЬхАТ, в) = 0.

Бх

(15)

^=М“а+М“шг + 8в,

(16)

Отклонение руля высоты задается в виде

5, при 1811 < 8,

8° 8maxSІgn81 при |81|>8тах,

8| - К, ^2'

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

8щах =: О И ^тах:=

В соответствии с моделью атмосферной турбулентности [7] величина является случайной для различных реализаций с плотностью вероятностей

о 2 »2

чг

2*2

где Р0= 1——Р2 — вероятность полета без ветра, Ри Р%, Ь1{ Ь% — заданные параметры.

Пусть найдена зависимость /?0(а. °(г0) Для фиксированной величины среднего квадратического отклонения скорости ветра о^, , тогда ДЛЯ любой Другой величины а\г очевидно будем иметь

аХГ

Аналогичным соотношением определяется и среднее квадратическое отклонение угла атаки

5,

.Ы«о.„(о,г0)-7-. (18)

О,

И7о

Если при фиксированной величине о-яг определена вероятность

п/ I , 1 —^о(®п7)Т’

Р(а | а иг) = 1 — е ^

(см. (7)], то полная вероятность с учетом случайности величины Оу,' определяется по формуле

00

Ре (ос) у {аур) Р \ йоур, (19)

о

Величина Рї представлена в виде

1-е V (20)

Вычисление значения Хе проводилось для интервала времени Т— = 1 час.

Таким образом, решение задачи сводится к следующей последовательности вычислений:

1. Определяется реакция системы ц4(0=а(0 на единичный импульс 5(0 =6(0), при фиксированной величине <п^= з^0.

2. Определяются величины за и та по формулам (12) и (13).

6„(0 = я

и проводятся расчеты для нескольких значений Я. Таким образом определяются зависимость ^0(а) и производные да/дЯ0.

4. Для фиксированного значения а* и соответственно Яо(а*) проводятся расчеты по методу Монте-Карло с формированием случайных реализаций |(0 по формуле (14) определяется величина йа — = Ж[(а— а*)2] и по формуле (15) находится значение функционала Я(а*).

5. По формуле (9) определяется величина %{оуг) и по формуле (20) с учетом соотношений (7), (17) и (18), (19) величина

Для вычисления величины 'Ао по формуле (9) использовались асимптотические выражения

У2

ЯЦа)

г 2 сіо. -.

тт.. Яо<**)

V 2*т„

1 ¿а

Х(о\г) =

2 а1

У 2 * т.

№

— (а*)

«Ю0

величины

■’«0»

'•«О

и зависимость Я0(а) определяются при фиксиро-

ванной величине

ОС. 10° м/с - 20

0 - 0

-10° - -20

20° - -40 ]

10° - 20

а - 0

-10° - -20

-20° - -чо

10° - 20

о - 0

Г - -20

20° - -40

град

10

Рис. 2

Рис. 1

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

= —0,4-|-; ма=-\^; УИ“ = 0;

Ь = 760 м, 1/ = 230 - , Кш = — 1 - •

’ с с

Величины Р1г Рг, Ьь Ь2 заданы для высоты Н= 11 км:

Р1== 0,01; />3 = 0,0001;

' ~ 1 м/с; Ь2^3 м/с.

Определение зависимости Яо(а) проводилось при значении а^о = 1м/с_

На рис. 1 приведены экстремальные реализации ветра №у(1) и со-

ответствующие им зависимости угла атаки а(£) от времени для /? = 56 и различных значений 6Шах-

На рис. 2 приведены зависимости Хе от величины а* для различных значений бшах- Штриховой линией показана зависимость Х^а*), определенная без поправки (15).

Отметим, что для приведенного примера (6тах= 1 град/с2) поправки (15) оказывают наибольшее влияние. Для предельных случаев

бтах = 0 И бтах =

+ оо поправки равны нулю {йа = 0).

ЛИТЕРАТУРА

1. Тихонов В. И. Выбросы случайных процессов. — М.: Наука,

1970.

2. Стратонович Р. Л. Избранные вопросы теории флуктуаций в радиотехнике. — М.: Советское радио, 1961.

3. К а з а к о в И. Е., Мальчиков С. В. Синтез стохастических систем в пространстве состояний. — М.: Наука, 1983.

4. Пугачев В. С., Синицын И. Н. Стохастические дифференциальные системы. — М.: Наука, 1985.

5. Вентцель А. Д., Фрейдлин Н. И. Флуктуации в динамических системах под действием малых'случайных возмущений. — М.: Наука, 1979.

6. Кузьмин В. П., Я р о ш е в с к и й В. А. Оценка предельных отклонений параметров траектории самолета при автоматической посадке.— Ученые записки ЦАГИ, 1984, т. 15, № 2.

7. Динамическое нагружение самолета при полете в неспокойной атмосфере. — Обзор ОНТИ ЦАГИ № 663, 1963.

Рукопись поступила 9/XI! 1988 г.