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

Том XXVI

1995

№3-4

УДК 629. 735. 33. 015. 073 519.2

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

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

Рассмотрим такую линейную стационарную систему:

где х — и-мерный вектор фазовых координат, % — /и-мерный вектор некоррелированных гауссовских белых шумов единичной интенсивности, А — матрица размера и х и, В — матрица размера и х т.

Если в начальный момент времени t — 0 вектор х имеет гауссовское распределение вероятностей, то он будет иметь гауссовское рас-

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

Приведена методика приближенной оценки вероятности первого пересечения гауссовским случайным процессом х (/) высокого уровня х.. Искомая вероятность определяется формулой Р (/) = 1 — е ^. Методика основана на определении среднего числа всех пересечений , среднего числа пересечений в одной серии п„ и представлении величины щ>

х = Ах + В£,

(1)

пределение при любом значении /, причем его среднее х (/) и корреляционная матрица К (Г) определяются соответственно уравнениями:

х = Ах, (2)

К = АК +КАТ + ВВТ. (3)

Пусть уравнения (2) и (3) имеют при / -> оо стационарные решения х = 0 и . Необходимо определить вероятность Р (/) того, что координата х\ (/) пересечет за время / уровень х. = стхй • Л, где ст1й — стационарное значение среднего квадратического отклонения координаты X]. Будем рассматривать случаи, когда величина Я достаточно велика, а начальное распределение вектора х таково, что Р (0) = 0. Тогда вероятность первого пересечения может быть представлена в виде [1]:

и основной задачей является определение величины

2. Дифференцируемый случайный процесс. Для упрощения последующих формул будем считать, что путем линейных преобразований переменных х-1, ... , хп система уравнений (1) приведена к виду

т

*1 =Х2 +2^7’

7=1

п т

*/• = £ аих] + £ ьу <4>

;=1 у-1

1 = 2,..., п.

Рассмотрим вначале пример, когда случайный процесс х\ (/) является дифференцируемым, т. е.

У-1

где ¿у — элементы первой строки матрицы В. В этом случае срденее ЧИСЛО пересечений ЭТИМ процессом уровня X. > 0 изнутри (Х1>0) в единицу времени определяется соотношением [1, 2]:

00

= \хр(х.,х,^с1х, (5)

о

где р (х, х, {) — совместная плотность распределения переменной *1 а ее производной.

Для стационарного гауссовского случайного процесса формула (5) примет вид

Разделим все реализации случайного процесса на две части: до момента первого пересечения и после пересечения. Очевидно, что второе и последующие можно рассматривать как пересечения, порожденные нестационарным случайным процессом, начинающимся в момент первого пересечения.

Первое пересечение характеризуется условием *1 = х. и условным распределением остальных фазовых координат р\ (х,, Х2, ..., х„) = = Р\ (хк)> которые являются начальными для нестационарного процесса, порождающего второе и последующие пересечения.

Пусть N (0 — среднее число пересечений уровня х» в единицу времени для исходного случайного процесса в момент времени /, а N1 (/, / + х) — среднее число пересечений в момент времени / + X для нестационарного процесса, начинающегося в момент времени и Тогда среднее число всех пересечений за время / определится соотношением

При больших значениях времени функция N± (t, t + т) зависит только от разности аргументов (t, t + х) = Ny (х), причем lim Ny (т) = 7Vst. Учитывая эти факты, перепишем соотношение (6) в

и продифференцируем обе части равенства по I. После формальных преобразований получим

При больших значениях времени / справедливы соотношения — — СІР

N (0« -Л^, —- = (¿о (1 - -Р(О)- В результате из формулы (7) следует,

jN (х)dx = j(х) 1 + jiVj (х, т + dx. (6)

о

о

т

виде

о

о

о

N (0 (1 - Р(0) = ^~(t) + (т) - ЛГЙ]*. (7)

о

что

о

00

где Л* = 1 +

о

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

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

Таким образов, задача об определении величины сводится к оценке функции N1 (х) и интеграла (8).

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

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

Р\ (Хк) — условное распределение^ фазовых координат в момент первого пересечения, а величина Ny(xR, х) определяется по формуле среднего числа пересечений для нестационарного гауссовского процесса ([2]): .

и, =

(*д)

(10)

Здесь

00

и»(*к) = 1 + ЦХ (хя, х) - ЛГЛ] ¿X

о

^і(хк,т) = ЛГ+(х) =

2п а-Дх)

(П)

где / (М+) = ехр - і МІ (х) + у/ЬІ М+ (х) Ф [М+ (х)],

где тх (х) и ту (т) — математические ожидания функции хх (т) и ¿1 (х) = у (х), которые определяются при решении уравнения (2) с начальными условиями х (0) = Хд; ах (х), ау (т), г (т) соответственно — средние квадратические отклонения переменной х\ и ее производной, и коэффициент корреляции между ними

и 2

Ф(и) = -,— Ге~в — функция Лапласа.

\2п

— 00

Корреляционная матрица для нестационарного процесса определяется при решении уравнения (3) с начальными условиями К (0) = = 0. Таким образом, из величин, входящих в формулу (11), от начальных условий х (0) = Хд зависят только средние значения тх (х) и ту (х).

Отметим также, что для гауссовского случайного процесса с математическим ожиданием х(т) и корреляционной матрицей К (/) известной является плотность распределения фазовых координат [6] р (х, /) и соответственно условная плотность распределения в момент пересечения уровня XI = х. для всех пересечений

р(ХяЛ=__^М1^п0У--------------------------. (12)

I ... 1р(х,,у,х3...х„,Оус1ус1х3...е1х„

— 00 0

Для нахождения распределения фазовых координат в момент первого пересечения (р 1 (хд)) используются следующие соображения. Нестационарный процесс, начинающийся от момента первого пересечения, порождает последующие пересечения — второе, третье и т. д., которые вместе с первым составляют все пересечения исходного случайного процесса. Поэтому функция р\ (хд) может быть определена из того условия, что она Ьместе с условным распределением фазовых координат в моменты пересечения уровня XI (() = х, для нестационарного процесса определяет известное распределение для всех пересечений исходного стационарного процесса.

Задача об отыскании распределения р\ (хд) в общем случае является весьма сложной и поэтому будет решаться приближенно. Распределение р\ (хд) используется для вычисления интеграла (10), поэтому потребная точность оценки этого распределения обусловлена зависимостью функции я, (хд) от начальных условий Хд. Если эта зависимость не очень существенна, то допустимы грубые оценки. В частности, если зависимость функции пх (хд) от переменной хд в окрестности среднего значения Хд близка к линейной, то соотношение (9) можно заменить на приближенное

л» = ¡п, (хл)р1 (хк)с1хя * л* (хк)

и ограничиться определением среднего значения Хд. Для этого используем соотношение, связывающее известное среднее значение

дсдй Для всех пересечений и средние значения хл и Хл„л соответственно для первого пересечения и для нестационарного процесса. Это соотношение имеет вид

где р (рс., у, /) — плотность распределения коородинаты Х\ и ее производной Х!=у, а ЗсЛ„51 — условное математическое ожидание вектора Хд при заданных значениях *1 = х, и = у для нестационарного процесса с фиксированными начальными условиями х (0) = Хд. Вектор Хд„&1 определяется соотношением [6]:

где компоненты векторов «1 (0 И «2 (0 определяются через элементы корреляционной матрицы К (() нестационарного процесса ([6]):

В соотношениях (14) и (15) от неизвестного вектора Хд зависит только вектор х (/), первые две компоненты которого обозначены как /Их(/) и ту{1). Соотношение (13) является системой п — 1 нелинейных уравнений для неизвестных хк2х^.

Интегрирование по переменной у в соотношении (13) с учетом соотношений (14) и (15) сводится к вычислению функций

(13)

1 + Цр(х. ,у, /)у4уЛ

о о

= X (0 + ах (0 (х, - тх (0) + а2 (0 (у - ту (¡)), (14)

а2г = (*/2 ■ *11 ~ *12 • */1) / (*22 ' *11 " *?2)>

аИ = (*/1 ~ *12 аИ ) / */1>

/ =1,2,..., Л.

(15)

00

о

00

о

о2уШ-г2«)) Г 1(Х9

—-------------ехр-------

2пах(() 21,

X

(16)

X 42п(М+ +1)Ф(М+) + М+е 2

Входящие в соотношение (16) величины определены выше [см. (11)].

3. Аналитические оценки. Предложенный метод позволяет решать задачи численно и получать аналитические оценки для задач небольшой размерности. Рассмотрим пример системы второго порядка:

х + 2кх + х = 2 4к % (/),

(17)

где % (/) — белый шум единичной интенсивности.

Приближенные оценки величины цо для к « 1 приведены в работах [3—5]. В работах [4, 5] приведены результаты численных экспериментов для данной задачи.

Стационарное распределение вероятностей переменной х(/) яв-

характеризуется параметрами х = тх =0,

ляется гауссовским

X = ГПу = 0, СТХ = СТу

и

1.

__ 1 2 /

двух уровней х. = ± R. Тогда Nst = — e~R '2

п

Будем оценивать вероятность пересечения

Nl(x) = N+(x) + N_(x),

где величина N_ (т) определяется по формуле (11), если в последней изменить знаки величин тх (т) и ту (т) на противоположные.

В данном примере ДЛЯ вычисления величины Но необходимо определить распределение только одной переменной — производной Хд в момент первого пересечения pi (хл). Условное распределение модуля производной для всех пересечений является релеевским:

/КМ) = |*л|е

(18)

При численных расчетах неизвестное распределение р\ (хд) заменяется дискретным, а при аналитических оценках принимается релеевским с неизвестным параметром ст = аХя :

А (|*д|) = ^ехР

2с1

(19)

Аналитические оценки основаны на упрощении формулы (11) при малых (к «1) и больших (к »1) значениях параметра к.

В случае к » 1 упрощение связано с тем, что уравнение (17) имеет два существенно различных корня характеристического уравнения: Xi и — 2к и Я.2 « — 1/2к. В результате изменение величин т^х),

(1\ ту(х), ах(х), ау(т) и г(т) при т > т, = О — J определяется прибли-

женными соотношениями:

тх (х) *

R +

XR 2 к )

2 к

2 к

о2х (т) * 1 - е Ф, с2 (т) * 1, г(т) * 0.

Величина п. при больших к и К определяется интегралом (8), который накапливается в основном на интервале 0 < х < х.., где х.. = = О(к/Л2). Поэтому если х.. » т. (или к2 » Л2), то формулы (20) можно распространить на весь интервал 0 < х < «. Кроме того, результаты расчетов зависимости п. (хц) показывают, что в этом случае распределение р\(хц) можно считать релеевским с параметром ст = 1.

Учитывая принятые допущения, из соотношений (8), (11) можно получить приближенное выражение для величины п, (хд):

ОО

exp

R

2к)

2(1 - e~%tk)

-exp

dx. (21)

Сравнение результатов расчетов, полученных по приближенной формуле (21) и по точным формулам (8), (11), показывает, что для х, > 0 они оказываются более точными при х < 0, чем при х > 0. Величина п, при х < 0 определяет [см. (9)] среднее время первого пересечения заданного уровня х, = R извне, которое при больших значениях R близко к среднему времени пересечения изнутри (9). В этом случае из соотношения (21) следует, что при больших к/Я

На.

V 7г R

Rx

R

2 к

Учитывая, что распределение производной xR принимается релеевским (19) с ст = 1, получаем:

о ГГ . г - Г- -Л

П*= I Pi (xR)n. (xR)dxR *l + JijeXp

(22)

При к « 1 упрощение соотношения (11) связано с тем, что функция тх (х) является колебательной и слабозатухающей. Очевидно, что пересечения уровней х. = ± R нестационарным процессом, происходят в окрестностях точек х = In, I = 0, 1, 2, ... . В этом случае изменение параметров нестационароного процесса при не очень малых х приближенно описывается формулами:

тх (х) = е~кх (i?cosx + xR sinx),

" W = °v (x)»l-e-2*V(x)*0.

Функции N+ (т) и N_ (т) имеют локальные максимумы в окрестностях точек х = Ы, когда среднее значение тх (т) оказывается близким к R или -R.

С учетом сделанных предположений найдем, что

М*д)-1 + £

/=1

Ф

R - ^R2 +x2r е~ы

4\-е^ы

-Ф (К)

(23)

00

_ 1 Г 2

где Ф(.К) = -= \е~и I1 <1и — функция, дополнительная к функции \2п ^

X

Лапласа.

При малых к сумму (23) можно заменить интегралом

со

М*д) = 1 + -[

71 J

Ф

R - jR2 +х2 е~кх

Ф (R)

dx.

Выполнив замену переменных и = 1-е т*, преобразуем этот интеграл к виду

И* (Хд) * 1 + J Ф

♦2 ^ « Х9

Ru - -4 2 R

• 2 ^

'-—2 4 R2

•Jlu

du.

При больших R значение п, (хК) определяется в основном окрестностью и = О, что позволяет получить оценку

и»(*д)*1 +

1 *R '

1+т

1-

4 R2

п kRi

(24)

Для случая к « 1 существенным оказывается отличие распределения производной (р\ (хд)) в момент первого пересечения от аналогичного распределения для всех пересечений (18). Определение параметра ст в распределении (19) основано на вычислении функции распределения производной х для всех пересечений.

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

где п. (у, хк) — среднее число пересечений для нестационарного процесса с начальным условием х (0) = хя и с производной в момент пересечения уровней ± Я, по модулю превышающей величину у.

Тогда безусловная функция распределения производной для всех пересечений определяется соотношением

Соотношение (25) является интегральным уравнением, которое можно решать численно, заменяя неизвестное распределение р\ (хЛ) дискретным и сводя уравнение (25) к системе линейных алгебраических уравнений.

Для грубых аналитических оценок распределение р\ (хк) примем в виде (19) с неизвестным параметром ст. Соотношение (25) при малых значениях кЛ и 1/Л2 можно преобразовать к виду

Еу (и) — модифицированная интегральная показательная функция,

оо

|Р(|л:| > у|хд|) рх {хК)йхК = е у2/2

(25)

О

Сокращая обе части равенства на общий множетель е у ^ и заменяя у2 на среднее значение у2, получаем

ст2 (1 + 1/Л2)

еи-Ех{и)

(26)

где

и = (1-ст2)(1 + 71М2)/ст2,

которая при больших и определяется разложением

Е\ («) =

е

(27)

Из соотношений (27) и (26) следует, что

(28)

Учитывая, что полученная формула справедлива лишь при малых значениях 1/Л2 и пкК2 и что при больших значениях пкЛ2 величина а стремится к единице, заменяем эту формулу эмпирической:

а2 =

1/Л2 + пкИ2

1 + пкВ2

(29)

На рис. 1 приведено сравнение значений а2, вычисленных по формуле (29) (штриховые линии), с результатами численных расчетов, в которых распределение р\ (х^) заменялось дискретным по 10 точкам (сплошные линии).

вл

0,8

0,7

0,6

цз

ол

10

Рис. 1

В итоге, используя формулы (24) и (29), получаем при малых к.

ґ г-:--------------\

00

«*=[«* (хЯ)Р! (хя)с1хл * 1 +-------------------;

ТГІГІ

пкЯ

1 + -

Л+пкк2

Л?______

1 + пкВ2

(30)

Для больших Л можно пренебречь членом 1/Л2 в (30) и объединить формулы (21) и (30) в одну:

пкЯ2

+ пкК2

(31)

Результаты расчетов по полученной формуле приведены на рис. 2. Сравнение полученных результатов с результатами числен-

ных расчетов из [4] показывает, что при Я = 3 погрешность формулы (30) не превосходит 20% при 1§ к й 0 (0,002 < к < 1) и 35% при ^ к > 0 (1 < к < 20). Можно предположить, что результаты численного эксперимента, полученные при к > 5, завышают значение п, из-за того, что при проведении такого эксперимента некоторые пересечения заданного уровня остаются незамеченными. Этого эффекта, по-видимому, можно избежать при использовании очень малого шага А интегрирования уравнения (17), удовлетворяющего при к > 1 усло-

, 1 вию А «—г-

2 к

При к » 1 или к «1 пересечения происходят сериями, и величина и, определяет среднее число пересечений в одной серии, следующей за первым пересечением. Приведенные на рис. 2 результаты показывают, что при фиксированном Л число пересечений в одной серии неограниченно растет как при к -> 0 так и при к -> <х>. С другой стороны, при любом фиксированном значении к величина п. стремится к единице при К -» оо.

4. Недифференцируемый случайный процесс. Рассмотрим теперь случай недифференцируемой переменной *!(/), при котором для сит

стемы уравнений (4) должно выполняться условие Ь2 = Ь}1 * 0.

7=1 1

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

. 1

г = Т

т

ЕМ/-*

и=1 .

и перепишем систему уравнений (4) в виде

*1 =х2 +Z = У,

г = ■

т

л

-г

и-1 .

т

X, = х} + ^Ьу = 2,...,и.

7-і 7-і

(32)

Переменная XI станет дифференцируемой, но при Тг -> 0 решение задачи для системы (32) будет стремиться к решению для исходной системы (1) или (4). Легко видеть, что

где / — 2, ... , п, гХ/Х1 — коэффициенты корреляции между переменными X/ и XI.

Представим соотношение (8) в виде

Формально величину (т) можно вычислять по формуле (11), если в последней положить ау = 1, гпу = 0 и г = 0. Функции тх{т) и стХ1(т), так же как и в случае дифференцируемой переменной, определяются в результате решения систем уравнений (2) и (3).

Величина й,, в соотношении (33) является размерной и имеет размерность производной Ж/сЬс. Поэтому в случае недифференцируемой переменной х (/) удобно ввести безразмерную величину й», определяемую по формуле

Величина й,, так же как и величина п. в случае дифференцируемой переменной, удовлетворяет предельному соотношению lim й, = 1. Действительно, при больших значениях R интеграл

<**! гх.хх -> 0 при Tz -> 0,

(33).

оо

Тогда, учитывая пределы (33), при Tz -» 0, получим

(34)

где

00

(35)

о

(36)

R-* 00

- 1 г

”* = 2я 1

■ “ 2' '

1 1 х, - тх (т) 1 е-*2/2

еХ1 (.V 1« 1 1 1 ст*1<т) ) С

ек

(37)

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

СТх, (Т) * ЬЪ> тх (т) * ту (°) • т’

*1 /иу(0)

^(0)

ьх

-Л.

2ст,

Подставляя эти соотношения в (37), получаем:

(38)

- 1 Г 1

Я**ТТ ТехР • 2пЬ1 * ут

1 Ь~ т 1 1 |>2

8 сР

Я* А

-я

ЯЬ:

Отсюда с учетом (36) следует, что й* « 1.

Рассмотрим пример системы уравнений второго порядка:

Х = у + е4ь 1

у = -2ку - х + 2-Лс

(39)

При е = 0 уравнения (39) сводятся к рассмотренному ранее уравнению (17).

В рассмотренном примере для вычисления параметра цо необходимо определить распределение переменной у в момент пересечения случайным процессом х (/) уровня х, = Я. При проведении расчетов величина й, определяется по формуле (35) для одного начального условия у(0) = уд, а условное среднее определяется по соотношениям, которые получаются из (13), (14) и (15^ предельным переходом при ст^ -> оо. В частности, в

этом случае в соотношении (15) компоненты вектора <*2 равны нулю. Величина й* определяется по формуле (36).

Результаты численных расчетов для различных зйачений параметров к, Я и г приведены на рис. 3, откуда видно, что величина й» значительно отличается от своего предельного (при Я -> оо) значения при малых значениях е. Для данного примера легко оце-Рис. з нить, что соотношения (38), при-

к=*5,0

нятые при оценке предельного значения й, = 1, выполняются, если Л » 1/е2.

5. Оценки вероятности выбросов угла атаки самолета. Рассмотрим короткопериодическое движение самолета под действием вертикальных порывов ветра ((). Линеаризованные уравнения возмущенного движения самолета вместе с уравнениями для формирующего фильтра (модель Драйдена) задаются в виде ([7])

с/а

1

СІСО

л

1 л

И?

ац,%.

(40)

В данном случае через а обозначена величина приращения воздушного угла атаки по сравнению с углом атаки, соответствующим

горизонтальному полету. Эта величина является недифференцируемой переменной.

Оценивается вероятность

Р (а (/) > а,), где а, = • Л. Для

больших значений Л величина вероятности определяется соотношением Р(а(0 > а*) « 1 - е~^а*.

Определение неизвестного распределения переменных юг, Ц'у и Z в момент первого пересечения заданного уровня угла атаки заменяется определением условных средних значений этих переменных.

При проведении расчетов принимались следующие численные значения параметров, входящих в уравнения (40): Ь = 760 м,

V = 250 м/с, М* = 1/з М“г. Варьировались частота О и относительное демпфирование к короткопериодического движения самолета и величина л“. При этом коэффициенты

М“ и определялись из соот-

Рис. 4 ношений

На рис. 4 приведены зависимости величины й* от Я, вычисленные по формулам (35), (36), для различных значений параметров и“,

П и к. Как и следовало ожидать, при больших значениях Я величина й, стремится к единице. Для рассмотренных значений параметров заметное отличие величины й* от единицы наблюдается лишь при малом относительном демпфировании (к < 0,2) и Я < 3.

Если принять, что й, = 1; то величина определяемая по соотношениям (36) и (34), составит

2 Lyfli

Vaaaj

^ 2 Re-*1'2.

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

ЛИТЕРАТУРА

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

2. Т и х о и о в В. И. Нелинейные преобразования случайных процессов,— М.: Радио и связь,— 1986.

3. Диментберг М. Ф. Нелинейные стохастические задачи механических колебаний.— М.: Наука.— 1980.

4. Кузьмин В. П., Ярошевский В. А. Оценка вероятности больших выбросов в динамической системе второго порядка при наличии случайных возмущений // Известия РАН. Техническая кибернетика.— 1993,

N»2.

5. Crandall S. Н. First-crossing probabilities of the linear oscillator //

J. of Sound and Uibration.— 1970, № 12 (3).

6. Вентце ль E. С. Теория вероятностей.— М.: Наука.— 1958.

7. Виноградов Ю. А., Ярошевский В. А. Выбор формы представления аэродинамических характеристик демпфирования летательного аппарата // Ученые записки ЦАГИ.— 1995. Т. 26, № 1 — 2.

Рукопись поступила 2/VIII1994 г.