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

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

№ 3

УДК 629.735.33.015.073

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

В. П. Кузьмин, В. А. Ярошевский

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

1. Постановка задачи. Рассмотрим динамическую систему, возбуждаемую случайным процессом белого шума на фиксированном интервале времени [*<,, Т]

xr=f{x, t)Jгg(x, *)Е, . (1)

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

Рассматривается случайный процесс х(0 с заданными начальными условиями х (/о) =Хо- Стохастическая система (1) понимается в смысле Стратовосича (1]. Уравнения подобного типа часто используются для решения задач динамики полета. Так, например, при анализе влияния атмосферной турбулентности на угловое и траекторное движение самолета можно аппроксимировать спектральную плоскость ветровых порывов дробно-рациональной функцией частоты (модель Драйдена, {2]). В результате к уравнениям движения самолета добавляются уравнения формирующих фильтров, в правых частях которых находится белый шум, а вектор фазовых координат расширяется.

Отличительная черта многих проблем динамики полета заключается в том, что выполнение поставленной перед летательным аппаратом задачи должно обеспечиваться с очень высокой степенью надежности. К таким проблемам относится, например, автоматическая посадка: согласно заданным нормативам вероятность летных происшествий не должна превышать Р% = Ю-6^-10~8 [3].

Это обстоятельство порождает новые подходы к решению задач, которые позволили бы найти гарантированные оценки для «выходных» параметров динамической системы [4, 5].

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

Будем считать выходным параметром случайное конечное значение координаты х1 (Т) и поставим задачу об определении вероятности Р* превышения некоторого заданного значения х*

Р [*1 (Г)>х,1 = Р*.

Если разбить рассматриваемый интервал времени [^о, Т] на отрезки с шагом At и заменить случайные функции кусочно-постоянными, то задача сведется к конечномерной. Оценки вероятности для динамической системы будут получены предельным переходом при Д£->-0.

2. Оценки вероятности для конечномерного случая. Рассмотрим исходные соотношения для конечномерного случая, когда возмущениями являются к независимых гауссовских случайных величин Си-.., с* имеющих единичную дисперсию и нулевое среднее. Пусть задана некоторая функция Хт=Хт {Си . . ., Си) и необходимо определить вероятность Р*=Р(*г>**)-

При малых значениях Р* приближенную оценку можно получить в результате решения задачи на условный экстремум [6, 7]

к

/?2 = тт V с? (2)

...°к ,=1

при условии хт (сь . . . , Си) = **• Оценка вероятности в первом приближении определяется соотношением

, 00

Рг = Р(ХГ>Х*) = ф (Я) = у=~Л. (3)

Если вектор с* является решением задачи (2), то формула (3) соответствует аппроксимации поверхности уровня Хт{с........ си) = х% каса-

тельной гиперплоскостью в окрестности точки с*.

Выполним ортонормированное преобразование (поворот) с=ВУ таким образом, чтобы вектору с* соответствовал вектор V* = = {#, 0, ... , 0}. Предположим, что поверхность уровня хт(У) = х% является гладкой, тогда в окрестности точки V* она может быть представлена в виде

у.=л + -эГ25>«1',и/. (4>

1=21=2

Очевидно, что касательная гиперплоскость определяется соотношением У1 = Р, что и дает первое приближение для оценки вероятности (3).

Ограничимся учетом кривизны поверхности уровня в форме (4) и будем считать, что второе приближение сводится к уточнению вероятности Р{хт>хц.) путем определения многомерного интеграла

+ СО +0С 03 1 2 ?

Р(*г>*,)=/ ■ • ■ /

—00 —00 у*

где

ъ=х+ж'к'к,,‘1у‘''>-

1—2 1=2

Предполагая, что Я достаточно велико и используя метод перевала для интегралов Лапласа [8], можно получить второе приближение для оценки вероятности

Р2 як Р(хт > х*)

(Я,)1'2

• С ~~7=г

1

(А)

1/2

Ри

(5>

где £>1 = с1е1:[/+Л], I — единичная матрица, А = [0,5] — матрица размерности п—1, составленная из коэффициентов, входящих в выражение (4).

Рассмотрим связь между матрицей А и матрицей, определяющей вторую вариацию функционала для задачи (2).

В соответствии с (4) представим функцию х(Уи...,Ук) в окрестности точки О,..., 0} в виде

х

дх

'* 1 дУх

(Уг -/?) +

1=1/=1

Используя множитель Лагранжа V запишем расширенный функционал для задачи (2) в виде

дх

у-Е1/‘ + -те

1=1

из условия 6/ = 0 получим

1=1/=1

Ж

дх

~дУ[

Тогда выражение для второй вариации примет вид

п к к

827= ^ 8 V? + £ ЕЬУ1 = 8У-,=1 /=1

(6>

;=1

Как следует из соотношения (6), матрица 1+А, используемая для оценки вероятности (5), совпадает с алгебраическим дополнением к первому диагональному элементу аи матрицы А, определяющей вторую вариацию функционала I.

Для дальнейшего существенным оказывается следующее правило. Если получено выражение для второй вариации функционала 1 от исходных переменных С\,..., Си\ ЬЧ = ЬсЧЬс и экстремальный вектор с*,. то для вычисления определителя в (5) необходимо выполнить упомянутый выше поворот системы координат с = ВУ и вычислить определитель алгебраического дополнения к первому диагональному элементу матрицы ВТ/ГВ.

При этом справедлива формула

(7>

где Д0 = (1е1|/:'|, а >.= тт 82/ при условии (8сс*) = /?.

...ск

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

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

Разобъем интервал [?о, Т] на N отрезков с шагом At. Функцию на /-м интервале представим в виде

«г

УМ

тп N

При Д*->0 условие = перейдет в условие

/=1/=1

т Т

5? (О л.= /?2

1=1 <0

и задача оценки вероятности Р{х1{Т) >х*) в первом приближении сведется к вариационной [6], а во втором приближении — к исследованию вариаций относительно оптимального решения («соседние экстремали» [9]).

Расширенный функционал для вариационной задачи имеет вид _ т Т

■/=Е[ + »(*1<7’)-**).

1=1 к

Экстремальные функции |г-(0 находятся из решения системы уравнений

х — г(х, t) + g(x, *)£=/(*, ?, 0.

-• _ дН дх

с граничными условиями

х(/0)е=х0; —

М7’) = *; \(Т) = 0, / = 2, ... ,я,

П

где Н= ^ Хг/г — функция Гамильтона.

/=1

Оптимальные значения переменных £г(0 вычисляются из условия

Е* (^) = агё Ш1П//(х, 5, *)•

Пусть найдено решение вариационной задачи £*(0, **(0> ^*(0> тогда вторая вариация функционала, обусловленная малыми вариациями б|(/), записывается в виде [9]

т

где

матрицы Н55 и /е определяются соотношениями

д*Н

дб/дбу

, /е =

ал

Матрица 5 удовлетворяет уравнению

5 = - 5/, —Гх 5 - Я,, + (Нл + 5/5) Щ1 (Я£, +/15) (8)

с граничным условием 5(7) =0.

Малые отклонения бх(^) описываются линеаризованным в окрестности экстремали уравнением (1) §х = Ух&х /е6£ и, следовательно, могут быть представлены в виде

I

Ъх(()= I* [!.(*, х)8|(т)йх,

4

где ц(£, т)—матричная импульсная переходная функция. Учитывая, что для данной задачи Я^ = 2/, запишем выражение для второй вариации в виде

§2

йі. (9)

Запись || а ||2 означает сумму квадратов элементов вектора [|а||2 =

Рассмотрим вначале случай скалярной переменной |(0- Введем обозначение

А (*, X) = 1 (*) + /1 (х) 5 (х)) ц (*, х).

Разобьем интервал [^о, Т] на отрезки с шагом At и обозначим через ij середину /-го отрезка. Вариации функции §(0 в дискретном виде за-

Ь с /

даются формулой 8^ = ^=-, а функция А (£, т) заменяется матрицей

^з=А У и ^).

Выражение (9) может быть представлено в виде

N

ЬЧ=^у], / = 1

где

У1 — ~2 -^и ^ + 8 С,;

Уъ = АпЪс1 М + Л22 8с2 М + 8с2.

Векторы у я 8с связаны соотношениями у = ВЬс, где В — треугольная матрица с диагональными элементами

Вп=тАи Д/+1.

Тогда 83У=8ст£тБ8с и определитель D() матрицы ВтВ определяется соотношением

В пределе при Д/-*0 получим О0 = ехр | АЦ, *) учитывая, что

*0

£) = /$(£), окончательно получим

£)0 = ехр

Определитель £>! вычислим по аналогии с формулой (7), которая при \t-y- 0 примет вид

П — —

^1- г ,

т

где Г=т1п§2/, при условии 18Ц*£^ = /?.

*5(0 4

10

Учитывая выражение для первой вариации функционала ЪJ= — ■h1(T)Ъx1 (Г) = 2{86(*) 5* (0^ = 2/?,

получим

XV -----------

МГ) •

8л:, (Л (Ю)

Определение величины А, сводится к решению задачи оптимизации квадратичного функционала (9) при линейном ограничении (10). Методы решения таких задач подробно описаны в [9]. Величина к определяется соотношением

Г=- 4 8*! (7Т в”1 (*0)8*! (Л, (11)

а функция 0(^) определяется в результате решения системы уравнений М = —/ТХМ + (Я,£ + 5/е) Нй1А М-,

Ь = М'/,Нй1АМ,

где М (£) — вектор-функция размерности т, а 0 — скалярная функция, удовлетворяющие граничным условиям

М1(Т)= 1, М,(Т) = 0, 1 = 2 6(Г) = 0.

Матрица 5(£) определяется уравнением (8). Оценка вероятности с учетом вторых вариаций определяется соотношением (5).

Если возмущением является вектор некоррелированных гауссовских белых шумов |(0> то по аналогии с предыдущим можно показать что

т

• Я0 = ехр [ [ 77? (1 (//«,+/* 5)/,= <«)],

%

где 77?(Л) — след матрицы Л.

Величина £>1( как и в случае скалярной функции |(0, определяется соотношениями (7), (11).

Рассмотрим теперь случай, когда кроме вектора белых шумов возмущениями являются и начальные условия для уравнения (1), которые определяются случайным вектором *о. Ограничимся случаем, когда вектор х0 имеет гауссовское распределение с нулевым средним и корреляционной матрицей КХо.

Представим вектор начальных условий в виде х0—Сс, где <ЗОт = Кх0, а вектор с имеет некоррелированные компоненты с единичной дисперсией.

Расширенный функционал для вариационной задачи имеет вид

п т Т

у=Е с?+ Е Г ??со ^(л-■*•).

1=1 /=1 у

*0

Решение задачи дает возможность найти вектор-функцию |*(0 и вектор с*. Вторая вариация функционала, обусловленная вариациями. б|(0 и 6 с, определяется формулой [9]

г

5* 5ст [I 5 (до + / ] 8с + 11 ж Ш ьг,

и

где

8гт = Н£ (Н^ + Д 5) 8х + 8$.

Определитель £>0 находится по формуле

т

Д, = бе1 От 5 (*0) О + /] ехр | +/| Я)/*) Л ].

Определитель вычисляется по формуле (7), где

X = пИп \ьх1 (Б — Ж0-1 + I) 8х0 + Ъх1 (Т) 6-1 М8х0 - 1 Ьх. (Т)2 —I

[ 2 ' 6 (*0) ]

при условии ХТ(/0)8л:0 —X! (Т)8д:1(7) = 2/?.

Введем п+ 1-мерный вектор Ъ\

гт={8*0^0), ЬХ1(Т)},

тогда определение величины X сводится к задаче квадратичного программирования с линейными связями X = minZтЛZ при условии Вл7-= 1, где

Г 5-7И6-1 М'+1 ЛГе-1'

А = I е-1^ - 6-1 ] *

£ = 27?1>^о), “МГ)]-

Все величины в матрице А берутся в момент времени /0- Решение задачи определяется соотношением

~В?~А-1В '

Оценка вероятности по-прежнему определяется формулой (5). 74

4. Пример. Рассмотрим применение полученных формул на примере одномерного уравнения

х = f(x) + У'21.

Функция f(x) имеет излом при |х\=а и задается формулой

j - х при \х\<а,

' 1—a sign х при |xj>a.

Рассматривается полубесконечный интервал [—со, Т] и граничное условие х(— оо) =0.

Необходимо определить вероятность Р(х(Т)^>х%).

Для данной задачи легко найти точное стационарное решение уравнения Фокера — Планка [9], которое позволяет определить плотность вероятности ю(х) и соответственно точное значение искомой вероятности

со

Я0(х(Г)>**)= [*{x)dx. (12)

Плотность вероятности очевидно удовлетворяет условию ш(л:) =со(—х) и для положительных значений х определяется формулой

»(■*) =

4-exp | — при 0<л: < а,

exp | — \{2ах — а3 + (1 + е) (х — a)2j

при х > а.

Нормированный множитель А определяется из условия

+ оо

2 j* a) (х) dx = 1.

о

Рассмотрим решение вариационной задачи для случая л:*>а>0.

Для полубесконечного интервала времени экстремали x*(t) удовлетворяют уравнению [61:

x=f(x). (13)

Оптимальная функция |(^) определяется соотношением (0= или

с учетом (13) S* (t) = ysif(x* (0).

Из приведенных соотношений легко получить

—ОС

В первом приближении вероятность Р(х(Т)>х%) определяется формулой

] — —

г“г■ <14>

я

Для получения оценки вероятности с учетом второй вариации необходимо определить функции S(t), M(t), Q(t), которые в данном примере удовлетворяют уравнениям

S = — 2fxS - 1ХЬ(х - a)s + S2;

M = —fxM + SB-,

0 = Л42

(15)

и граничным условиям

S(T)~> О, М(Т)= 1, 0(7’) = О,

где

J — 1 при |*|<а,

I — (1 4-е) при |х| >а,

б (л: — а) — дельта-функция.

Уравнения (15) решаются аналитически. Выберем систему отсчета времени таким образом, чтобы момент ^ = 0 совпадал с условием х*(0)=а. Тогда функция S(,t) претерпевает разрыв при t = 0:

О при t > 0,

s(0 = | 2S0e2t

21

при £<0,

1 — 50е

постоянная 50 определяется из условия

2 Зо —Ое

1 — 50 *

используя 5(^) можно последовательно получить формулы для опре-деления функций М(0 и 0(0-

Величина Д, определяется соотношением

г

, = exp J 5 {t) dt = 1 + s,

D —

где

2 R

2 6(— ooj ’ Ьх(Т)=* ЦТ) — a-ril + tHX' — a)

Выражение для 0(—оо) в общем случае оказывается очень грамозд-ким, поэтому рассмотрим некоторые частные случаи:

1) е>0. В этом случае при больших х% можно пренебречь изменением функции М(Ц 0(^) при /<0, так как

М (0) =

■**(! + е) —еа

и при л* >й М(0)<С1, а при ^<0 M(t) убывает быстрее, чем М(0)е Тогда

40):

А*, (Л:

1(—оо)~-

1

2(1 + *) ’

2 т1№'

При больших значениях х.м/а можно получить приближенные отношения оценок Р1 и Рг к точному значению Р0:

Р\ Ро

V 2г.

к 1 + е;

Результаты расчетов по формулам (14), (16), (19) приведены на рис. 1. Из приведенных результатов расчетов и приближенных формул видно, что в случае е>0 оценки Р1 и Р2 имеют одинаковый порядок точности.

2) е = — 1. В этом случае А>=0 и 1 = 0, поэтому величину следует определять как предел

г

= Нт ---р--------------.

<=“°°---5“ Лж» ( 7*)* в~‘ (<)

Для случая е = —1 легко получить

Ц,(0-|5 (<)<« = *";

*

в(0 = 4-(1-е-аО;

Ьх(Т) =

А =

_5_ а

а_\ 2 Я

Отношения оценок вероятности Р4 и Р2 К точному значению Ро) будут цри больших значениях х*/а определяться приближенными формулами

г-1,0

0,5

Рис. 2

Рис. 1

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

Отметим также, что при а> 1 нормировочный множитель Д ^ У2 те и оценка Р2 является практически точной.

ЛИТЕРАТУРА

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

2. Доброленский Ю. П. Динамика полета в неспокойной атмосфере.— М.: Машиностроение, 1969.

3. Белогородский С. Л. Автоматизация управления посадкой самолета. — М.: Транспорт, 1972.

4. Малышев В. В., К и б з у н А. И. Анализ и синтез высокоточного управления летательными аппаратами. — М.: Машиностроение,

1987.

5. К е й н В. М. Оптимизация систем управления по минимаксному критерию. — М.: Наука, 1985.

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

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

8. Федорюк М. В. Асимптотика. Интегралы и ряды. — М.: Наука, 1987.

9. Брайсон А., Хо Ю - Ш и. Прикладная теория оптимального управления. — М.: Мир, 1972.

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

Рукопись поступила 29/V 1989 г.