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

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

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

Работа выполнена при финансовой поддержке РФФИ, проект № 11-01-96024.

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

Затрудняет подобные исследования и то, что в отличие от стохастических дифференциальных уравнений (СДУ) с винеровскими возмущениями, для которых построение соответствующего дифференциального уравнения (в частных производных) Фоккера-Планка-Колмогорова (ФПК-уравнения) для плотности распределения фазового вектора системы (и для СДУ Ито, и для СДУ

Стратоновича) не составляет труда [3], добавление скачкообразного пуассоновского возмущения превращает ФПК-уравнение в интегро-дифференциальное уравнение (в частных производных) Колмогорова-Феллера (КФ-уравнения) [4], связь коэффициентов которого с системой СДУ не столь очевидна из-за различных форм пуассоновских процессов [5]. Не столь очевидна и форма стационарного решения КФ-уравнения даже в простейших случаях. Более того, даже для простейшей линейной системы с аддитивными смешанными возмущениями стационарная плотность вероятности неизвестна.

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

1. Постановка задачи

Рассмотрим систему, описываемую следующим стохастическим дифференциальным уравнением Ито:

dX ( t ) = f ( X ( t ) ,t ) dt + g ( X ( t ) ,t ) dW ( t )+

+ c ( X ( t ) ,t ) dP ( t ) , t>t о , X ( t о )= X о . (1)

Здесь t - время. X ( t ) - случайный процесс, характеризующий поведение исследуемой системы, X о - случайная величина с известной плотностью вероятности p о( x ), W ( t ) - стандартный нормальный винеровский случайный процесс [8]: E[ W ( t )] = 0, E[ W ²( t )] = t. P ( t ) - составной случайный пуссоновский процесс [8] (его характеристики уточним ниже), f ( •, • ), g ( •, • ), c ( •, • ) -неслучайные функции, E - символ математического ожидания.

Конкретизируем вид исследуемой системы, а именно предположим, что процесс P ( t ) представляет собой последовательность случайных импульсов прямоугольного вида постоянной ширины [1]:

N ( t )

P ( t ) = ^TZi • 6 ( t - ti ) ,          ( 2 )

i =1

где Zi - случайные амплитуды импульсов, имеюшие одинаковую плотность вероятности p* (z), причем импульсы появляются с постоянной интенсивностью X > 0. Кроме того, предположим, что f (x, t) = —a x, g(x, t) = в, c(x, t) = Y, где a,Y > 0, в - постоянные. В этом случае интегро-дифференциальное уравнение Колмогорова-Феллера [1,5] примет вид

dp ( x,t )     в ² д ² p ( x,t ) д [ axp ( x,t )]

dt      2 dx ²   + dx

—Xp ( x, t ) + X

+ ж

J p ( x — y z,t ) p* ( z ) dz =

-∞

= L[ p ( x,t )] , p ( x,t о )= p o( x ) , где p ( x, t ) - одноточечная плотность распределения случайного процесса X ( t ).

Ограничимся вычислением стационарной плотности ps ( x ). удовлетворяющей уравнению

L[ ps ( x )] = 0                (3)

и условиям

ps ( x ) > 0 ,

+с ю

У ps ( x ) dx = 1 ,

-∞

для случая нормально распределенных ам

плитуд импульсов:

p* ( z ) = N ( a*,^) =

2 Пк С*

exp

-

( z — a* )² 2 a ²

a* = 0 , (1)

где a* и a* >  0 - математическое ожидание и среднее квадратическое отклонение случайных величин Zi соответственно.

Исходя из имеющейся структуры уравнения (3) решение задачи будем искать в виде смеси нормальных гауссовых распределений, а неизвестные параметры этих плотностей вычислим на основе применения метода коллокаций (проекционного метода решения

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

2.    Схемы расчетов и их реализации

Перейдем в несобственном интеграле в уравнении (3) к переменной интегрирования q q = x - yz, z = x—q,    dz = -dl

γγ

Тогда уравнение (3) примет следующий вид:

в ² д ² ps ( x ) + д [ axps ( x )] 2 dx ²         dx

— ^Aps ^{( x ) +}

+ A

+ ∞

У Ps ⁽ q ) p* ⁽ q ) dq = o ,

-∞ где

p* ( q )= N ( x,a 2 ) , a * = ца*.     (6)

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

KK

"ps ( x )= 52 rk N ( ak ’ak) = 52 rk psk (x ), k=-K             k=-K где

K

52 rk = 1, rk ^ 0, r-k = rk, k=-K a 0 = 0, a-k = —ak, a- k = ak, k = 0,1, 2 ,...,K, что для полного определения аппроксимации плотности ps (x) приводит к необходимости выбора параметров rk, ak, k = 1, 2, ..., K. каким-либо методом ('значения ak считаем фиксированными, задавая их, например, как ak = k • h- h > 0 - постоянный шаг. ii.tii каким-либо иным способом).

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

K

Ф[ ps ( x )] = Е rk L[ psk ( x )] .      (7)

k = -K

Расчет же L[ psk ( x )] в данном случае может быть выполнен полностью, если воспользоваться соотношением для свертки двух нормальных распределений:

N ( m 1 , s 1) *N ( m 2 , s 2) = N ( m 1 + m 2 ,s 1 + s 2) .

С учетом этого получим, что ф = Е rk[{^[^ — > k=-K            σk      σk

+ _а ^[1 —          ] —A } N ( ak,ak ) +

+ AN ⁽ ak,ak ⁺ a 2)] . (8)

1°. Метод коллокаций. Согласно общей схеме этого метода необходимо выбрать некоторое множество узлов {xi} (например, как xi = h*i. h* > 0 - постоялиый шаг. i = 0. 2..... с учетом четности ps(x) в дан ной задаче), число которых совпадает с количеством параметров rk и ak, требующих определения, сформировать систему трансцендентных уравнений вида

Ф ^ ( xi )] =0             (9)

и решить ее относительно параметров.

Данная схема была реализована в среде пакета Mathematica, причем все необходимые выкладки по построению системы уравнений (9) были выполнены точно. Но попытка прямого использования функции FindRoot для приближенного решения построенной системы уравнений оказалась неудачной вследствие невозможности указания при вызове этой функции дополнительных ограничений на параметры:

0 <  rk < 1 , ak >  0 .         (10 )

Поэтому пришлось воспользоваться процедурой FindMinimum пакета Mathematica для поиска условного минимума функции

2 K +1

F ( г , а )= £ Ф2[ ps ( xi )]      (П)

i =0

с условиями (10). Заметим, что на этапе аналитических выкладок параметр r о был сразу же исключен из уравнений на основе учета нормировочного соотношения для плотности ps ( x ), что привело к уменьшению на единицу числа параметров, требующих определения.

На рис.1, 2 приведены некоторые результаты расчетов при следующих значениях констант (номера списков параметров соответствуют номерам кривых на графиках):

h = 1 / 8. h* = 1 / 16. K = 32. y = 1- A = 1 (для всех),

• 1)    a = 2. в = 2 / 10. с* = 6 / 8:

• 2)    a = 2. в = 2 / 10. а* = 12 / 8:

• 3)    a = 2. в = 2 / 10. а* = 20 / 8:

• 4)    a = 2. в = 1 / 8- а* = 20 / 8:

• 5)    a = 2. в = 1 / 5- а* = 20 / 8:

С) a = 3 / 4. в = 1 / 8. а* = 16 / 8.

2°. МНК. Общая схема этого метода основана на минимизации интеграла от квадрата невязки

H ( r , a

+ ∞

)=У Ф²[ p s ( x )] dx

∞

по векторным параметрам r и ст , причем используя соотношения

+ ∞

У N ( m 1 ,s 1) N ( m 2 ,s 2) dx =

-∞

\ Пк S 0

exp

[-

( m 2 - m 1)²

² s 0

,

vν

+ ∞

/ x ^N ( m,s ²⁾ dx.

v 1 = m.

-∞ v 2 = m2 + s2.  v з = m3 + 3 ms2.

v 4 = m ⁴ + 6 m ² s ² + 3 s ⁴ .  s 2 = s 1 + s 2 .

можно получить явную форму функции H(r. a) и попытаться найти ее минимум. К сожалению, даже при небольшом K (порядка 8.. 10) объем аналитических выкладок, время получения результата - вычисления минимума функции H (r. a) при условиях (10) (это время явным образом зависит от сложности H) и затраты оперативной памяти компьютера (для расчетов использовался ноутбук Samsung с процессором Intel Core Duo Т2450 2,0 ГГц 2,87 Гб ОЗУ) оказываются столь велики, что говорить об эффективности рассматриваемого метода для решения поставленной задачи не приходится. И это несмотря на то, что в схеме МНК отсутствует необходимость выбора множества узлов {xi}, что уменьшает неопределенность и, из общих соображений, должно повысить точность результатов.

Заключение

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