CC BY
47
УДК 519.63 + 519.246
ББК 22.161.6 + 22.193 + 22.171
СПЕКТРАЛЬНЫЙ МЕТОД АНАЛИЗА СТОХАСТИЧЕСКИХ СИСТЕМ С РАЗРЫВАМИ ТРАЕКТОРИЙ, ОПИСЫВАЕМЫМИ СЛУЧАЙНОЙ СМЕСЬЮ ЭРЛАНГОВСКИХ РАСПРЕДЕЛЕНИЙ1
Кожевников А. С.2, Рыбаков К. А.3
(Московский авиационный институт, Москва)
В статье рассматриваются стохастические системы управления с импульсными воздействиями, которые образуют гипер-эрланговские потоки событий и приводят к разрывам траекторий системы. Решается задача нахождения плотности вероятности вектора состояния. В основе решения лежит использование спектральной формы математического описания систем управления.
Ключевые слова: гиперэрланговское распределение, задача анализа, скачкообразный процесс, спектральный метод, стохастическая система.
1. Введение
Многие окружающие нас явления и закономерности (природные, технические, экономические) имеют случайный характер, что позволяет их описывать с помощью случайных процессов. Для описания моделей явлений, учитывающих влияние случайных факторов, как правило, применяют стохастические дифференциальные уравнения (при условии, конечно, что время можно считать непрерывным).
1 Работа выполнена при поддержке РФФИ (проект №12-08-00892-а).
2 Александр Сергеевич Кожевников, аспирант (exequit@yandex.ru).
3 Константин Александрович Рыбаков, кандидат физико-математических наук, доцент (rkoffice@mail.ru).
Используя уравнения с диффузионной и скачкообразной компонентами, можно моделировать поведение довольно сложных систем, учитывая как непрерывные случайные воздействия, так и импульсные, приводящие к разрывам траекторий. Однако часто для описания скачкообразной компоненты ограничиваются общим пуассоновским процессом (или пуассоновской случайной мерой) [19, 23, 24]. Такой выбор обусловлен развитым математическим аппаратом (стохастические дифференциальные уравнения с пуассоновской составляющей, уравнения Колмого-рова-Феллера), но следствием выбора общего пуассоновского процесса является использование только показательного закона распределения для промежутков времени между последовательными импульсными воздействиями (разрывами траекторий). В частности, для случайных величин, имеющих показательное распределение, среднее значение и стандартное отклонение совпадают, что вносит некоторые ограничения при решении прикладных задач [6].
Для того чтобы снять такие ограничения, можно использовать модели систем со случайным периодом квантования -частный случай систем со случайной структурой [4, 7, 12, 13, 20]. В них появлением разрывов траекторий управляет вспомогательный марковский процесс с конечным множеством состояний, задаваемый начальным распределением и интенсивностями переходов. Использование таких моделей приводит к ряду трудностей, в основном связанных со сложностью решения уравнений для нахождения плотности вероятности вектора состояния - системы обобщенных уравнений Фоккера-Планка-Колмогорова (интегро-дифференциальных уравнений в частных производных), решение уравнений для моментов тоже сложнее, чем для моделей, использующих только общий пуассоновский процесс. Вторая сложность связана с тем, что большое количество параметров модели (интенсивности, параметры распределений приращений, которые получает вектор состояния при разрывах траекторий для различных переходов между состояниями вспомогательного марковского процесса) затрудняет их идентификацию. Поэтому целесообразно, взяв за основу модель
системы со случайным периодом квантования, рассмотреть некоторые частные случаи с обозримым числом параметров. Один из таких частных случаев сводится к описанию промежутков времени между последовательными импульсными воздействиями эрланговским законом распределения, т.е. когда состояния вспомогательного марковского процесса последовательно сменяются «по кругу», а переход из некоторого, заранее выбранного, состояния в следующее сопровождается разрывом траектории: моменты появления разрывов траекторий образуют эрланговский поток событий. Такие модели изучены, например, в [4, 10, 22]. Модель, рассматриваемая в этой статье, допускает, что моменты появления разрывов траекторий образуют гипер-эрланговский поток событий и формируется два типа приращений при разрыве траектории, имеющих различные законы распределения. Распределение промежутков времени между последовательными разрывами можно задать смесью эрлангов-ских распределений. Использование такого гиперэрланговского распределения позволяет включить показательное и эрлангов-ское распределения как частные случаи и расширить множество решаемых прикладных задач. Аналогичный подход применен в [22, 25].
В последнее время системы, описываемые стохастическими дифференциальными уравнениями с диффузионной и скачкообразной компонентами, все чаще используются для моделирования цен акций [9, 21, 24]. Разумеется, задачи финансовой математики составляют только часть прикладных задач, при решении которых целесообразно использовать математический аппарат теории систем со случайным периодом квантования. Например, можно рассматривать асинхронные системы управления сложными техническими объектами [4], механические системы с импульсными воздействиями [22], изучать влияние случайных импульсных воздействий и помех на электрические цепи [2, 19] и др.
Все существующие методы анализа систем со случайной структурой, допускающих при смене структуры разрывы траекторий, можно применять и при анализе систем со случайным
периодом квантования. Один из них - метод Монте-Карло [1, 8]. Он сводится к многократному непосредственному моделированию реализаций случайного процесса и их статистической обработке. Существенное преимущество этого метода состоит в простоте реализации и универсальности. К его недостаткам можно отнести необходимость проведения большого количества вычислений для получения достаточной точности, особенно при оценке плотности вероятности вектора состояния.
Другие методы, так или иначе, связаны с уравнениями для нахождения плотности вероятности вектора состояния - обобщенными уравнениями Фоккера-Планка-Колмогорова: методы, основанные на задании структуры плотности вероятности вектора состояния, например, полигауссовская аппроксимация; переход к моментам; методы, связанные с представлением плотности вероятности или связанной с ней характеристики в виде функционального ряда: степенного, ортогонального; численные методы). Подробное описание таких методов содержится в монографиях [7, 20].
Предлагаемый метод решения основан на ортогональных разложениях функций с использованием спектральной формы математического описания систем управления [12, 18], он является развитием подхода, применяемого для более простых стохастических систем [2, 10, 11].
Отметим, что представление решения задачи в виде ортогонального ряда довольно часто используется при построении приближенно-аналитических методов. Но, как правило, для этого выбирается конкретная система ортогональных или биор-тогональных функций и выводятся соотношения для нахождения коэффициентов ряда. При таком подходе, конечно же, используются свойства выбранных функций (например, рекуррентные соотношения, связывающие эти функции; формулы, устанавливающие связи функций и их производных и т.п.). Основное отличие предлагаемого метода состоит в использовании произвольной ортонормированной системы, при этом соотношения для определения коэффициентов ряда не зависят от выбранных функций. Они представляются матричными уравне-
ниями, для которых можно получить формулы точного решения или сформировать методику приближенного решения. Свойства выбранных функций также используются, но несколько иначе, а именно от этих свойств зависит структура матриц, которые ставятся в соответствие линейным операторам (операторам умножения, дифференцирования, интегрирования и др.), определяющим исходную задачу. Выбор той или иной системы ортонормированных функций может привести к тому, что соответствующая матрица будет треугольной, ленточной, симметрической или кососимметрической. Для многих базисных систем, например, для полиномов Лежандра, тригонометрических функций, функций Уолша и Хаара, обобщенных функций Ла-герра (включающих полиномы и функции Лагерра) и обобщенных функций Эрмита (включающих полиномы и функции Эрмита), такие матрицы известны [12, 14, 16-18].
Для рассмотренной модели стохастической системы с разрывами траекторий, описываемыми случайной смесью эрлан-говских распределений, получено аналитическое решение матричных уравнений для коэффициентов разложения функций, по которым определяется плотность вероятности вектора состояния. Это же решение можно использовать и как приближенное, получая конечное число коэффициентов разложения и представляя плотность вероятности в виде частичной суммы ряда. Для апробации предложенной методики решен пример задачи анализа одномерной нелинейной стохастической системы с двумя типами приращений при разрыве траектории: случайным гауссовским и детерминированным; промежутки времени между последовательными разрывами описываются эрланговским и показательным распределениями.
2. Постановка задачи
Будем предполагать, что поведение модели системы управления описывается случайным процессом в непрерывном времени, являющимся аддитивной смесью диффузионного и специ-
ального скачкообразного процессов и удовлетворяющим стохастическому дифференциальному уравнению Ито [13]:
(1) dX(0 = f(t,X(t))dt + о^,Х(^^(^ + dQ(t), X(t0) = Х0, в котором X е Л” - вектор состояния; t е Т, Т = [t0, ^] - заданный отрезок времени функционирования системы;
У(/, х): Т х Л” ^ Л” - вектор-функция размеров ” х 1; о(^ х): Т х Л” ^ Л”х* - матричная функция размеров ” х 5;
Ж(0 - 5-мерный стандартный винеровский процесс, который не зависит от начального состояния Х0, определяемого заданной плотностью вероятности ф0(х).
Случайный процесс Я(() представляется в виде J а)
) = 1(6^,-) + (1 - 6 Шъ)),
г=1
где ./(0 - гиперэрланговский процесс, ассоциированный со случайным потоком событий [4], состоящих в том, что вектор состояния X получает приращения У\(т) е Л” или Г2(тг) е Л” в случайные моменты времени ть т2, ... Случайный вектор У\(т,) характеризуется плотностью вероятности q\(t, у), а случайный вектор У2(т,) - плотностью вероятности q2(t, у); t = тг. Выбор приращения У\(т) или Г2(тг) зависит от случайной величины Ъ>, принимающей значения 1 с вероятностью Х\/(Х\ + Х2) и 0 с вероятностью Л.2/(Л.1 + Х2) (случайные величины £1, £2, ... имеют распределение Бернулли):
X(тг) = X(тг - 0) + ^(Тг), 6 =1 }.
№(тг), 6 = 0}
Заданные положительные числа Х1 и Х2, а также натуральные числа N1 и N определяют гиперэрланговский закон распределения промежутков времени АтI = тг — гг-_1, г = 1, 2, ... (т0 = t0), который является эрланговским с параметрами Х1 и N1, если £ = 1, или эрланговским с параметрами Х2 и N, если £ = 0. Случайные величины £ независимы, поэтому выбор закона распределения для случайного приращения - ql(t, у) или q2(t, у) - в момент времени тг не зависит от предыстории.
Возможно и другое описание схемы появления событий (разрывов в траекториях вектора состояния). Для этого рассмот-52
рим случайный процесс К(^ с конечным множеством состояний {1, 2, ..., Щ, где N = N1 + Ы2 — 1. Интенсивности переходов задаются следующим образом (см. рис. 1): смена состояний 1^2, 2^3, ..., N1 — 1^М и N1^1 происходит с интенсивностью Х1, а смена состояний 1^-^ + 1, N1 + 1^-^ + 2, ..., N — 1^N и N ^1 - с интенсивностью Х2; другие переходы невозможны. При переходе из состояния с номером N1 в состояние с номером 1 вектор X получает случайное приращение Г1(тг), а при переходе из состояния с номером N в состояние с номером 1
- случайное приращение У2(т), что соответствует разрыву (скачку) траектории процесса X(t) (см. рис. 2).
Рис. 1. Граф состояний случайного процесса К(і)
Если предположить, что Х1 > 0, Х2 = 0, то /V) будет эрлан-говским процессом порядка N1, а процессу К(^ будет соответствовать кольцевой граф состояний (аналогично при Х1 = 0, Х2 > 0). Этот случай рассмотрен в работе [10].
Заметим, что постановка задачи может быть изменена с предположением, что интенсивности Х1 и Х2 зависят от времени, т.е. Х1 = Х1(/), Х2 = Х2(/). Далее будем рассматривать именно этот вариант.
Предполагается [24], что функции ^, х) и g(t, х) = а(), х) х х о-т(/, х) кусочно-непрерывны по t для всех X е Л”; при фиксированном t е Т функция /(1, х) имеет непрерывные и ограниченные частные производные первого порядка, а функция gij(t, х) -непрерывные и ограниченные частные производные второго порядка по координатам вектора состояния, г,у = 1, 2, ..., ”. Кроме того, Х1(/) и Х2(/) непрерывны на множестве Т, а случайные величины, задаваемые плотностями вероятности д^, у) и
Ч2(/, у), имеют конечные моменты второго порядка.
Рис. 2. Выборочные траектории процессов ^) и X(t)
Эти условия (иногда к ним добавляется дополнительное условие на матрицу диффузии g(t, x), требуемое для разрешимости приведенных ниже параболических уравнений) обеспечивают существование и единственность сильного решения уравнения
(1), однако они вносят слишком много ограничений при решении прикладных задач. Вопросы, связанные с ослаблением приведенных условий, рассмотрены в [3, 5]. В частности, можно рассматривать уравнения с разрывным коэффициентом сноса x) или вырожденной матрицей диффузии g(t, x), достаточно часто встречающиеся в задачах управления. Отметим, что в [3, 5] изучаются стохастические дифференциальные уравнения без скачкообразной компоненты или с пуассоновской составляющей, тем не менее эти результаты могут быть обобщены и на рассматриваемый случай. Кроме того, можно понимать решение уравнения (1) в слабом смысле, тем более что далее ставится задача нахождения плотности вероятности вектора состояния.
Таким образом, рассматривается стохастическая система с расширенным вектором состояния, непрерывная часть которого
- X, а дискретная - К є {1,2,...,Ж}. Тогда плотность вероятности ф(1, х) вектораXможет быть представлена в виде суммы:
N
р(1, х) = ЁРк -'О1, х),
к=1
где функции ф(к)(1,х) удовлетворяют системе обобщенных уравнений Фоккера-Планка-Колмогорова [4, 7, 12, 20]:
——(1, х) = Ар(1) (1, х) - ( Х (1) + Х (1)) р(1) (1, х) +
(2) 51
+Х(1) | q1(t,х - z)р(Д1)(1,z)dz + Х(1) | q2(t,х - z)р(д)(1,z)dz,
я" я"
5р—(1, х) = Ар(к) (1, х) - Х (1 )р(к) (1, х) + Х (1 )р(к-1) (1, х),
(3) 51
к = 2,..., Ду,
(4) 5р( 1 )(1, х) = Ар{ д+1 (1, х) - Х (1 )р( д+1 (1, х) + Х (1)р(1) (1, х),
51
5р—(1, х) = Ар(к) (1, х) - Х, (1 )р(к) (1, х) + Х, (1)р(к-1) (1, х),
(5) 51
к = Д1 + 2,..., Д, в которой
" 5 Г" “і
Ар(к) (1, х) = -Ё — Г у; (1, х)р(к) (1, х) ] +
;=1 5х; -1
(6) 2
і п п 52 г
+т Ё Ё ^(1, х)р(к) (1, x)], к = 1,2,..., д,
2;=1 у=1дх1дх] -1
gi]■ (1, х) = ЁСтг-Г (1,х)о]Г (1,х), і, у = 1,2,...,п.
г=1
Для случайного процесса К(1) начальное состояние фиксировано: К(10) = 1, поэтому с учетом заданной плотности вероятности ф0(х) начального состояния Х0 имеем
(7) р(1)(10,х) = р0(х), р(к)(10,х) = 0, к = 2,...,N.
Последние два слагаемых в правой части уравнения (2) для дальнейших преобразований целесообразно представить в операторной форме, определив
(8) Hrq)(t, x) = Xr (t) J qr (t, x - z)y(t, z)dz, r = 1,2,
Rn
для всех допустимых функций q(t, x); Hr - линейный оператор, а именно композиция оператора умножения на функцию Xr(t) и оператора Фредгольма с ядром qr(t, x - z).
Задача анализа стохастической системы управления, описываемой уравнением (1), заключается в нахождении плотности вероятности ^(t, x) вектора состояния X при заданных функциях ft, x), o(t, x), Ai(t), ^(t), qi(t, y), q2(t, y) и p0(x). Далее предполагается, что это решение существует в классе функций W21,1(T x Rn), т.е. плотность вероятности ^(t, x) и ее производные первого порядка принадлежат пространству L2(T х Rn), при этом решение {^(1)(t, x), ^(2)(t, x), ..., q>(N>(t, x)} системы (2)-(5) понимается в обобщенном смысле, т.е. эта совокупность удовлетворяет соотношениям
J ^(1)(t, x)——(^ x) dtdx = J <p(1)(t, x)A*y(1\t, x)dtdx -
T xRn ^t T xRn
- J i^(1)(t, x) (\(t) + ^2(t))^(1)(t, x)dtdx +
T xRn
f
+ J ^(1)(t,x) A1(t) J q1(t,x- z)q)('N1\t,z)dz +
T xRn V Rn
+ X2(t) J q2(t,x - z)^(N)(t,z)dz
Rn у
J y(k)(t, x) ——dtdx = J <p(k)(t, x)A*y(k)(t, x)dtdx -
T xRn ^t T xRn
J ^(k)(t,x)(^(t)^(k)(t,x)-^(t)^(k-1)(t,x))dtdx,
T xRn
k = 2,..., N1,
I х) <Д1х)*ь=
= | р<Д+1)<^, х)АУ<Д1+1)<^, х)&ёх -
Т xRn
- | у(Д1+1)<^,х)^(Ор1^1+1)<^,х)-Л^Ор1'1'1^,х))dtdx,
Т xRn
| у<к)<t, х) ——<t, х) dtdx = | р<к)<t, х)А*у<к)<t, x)dtdx -
Т xRn ^ Т xRn
- | У^О1, х)(л^Ор1^)<^ х) - Я^Ор^-1^, х))dtdx,
Т xRn
к = N. + 2,..., Д,
для любых бесконечно дифференцируемых функций ^'(1'l(t, х), ц(2:1<1, х), ..., Ц <д><1, х) с компактными носителями, где А* - сопряженный оператор:
Л><к % х) = 1/ <t, х) +
г=1 5гг
1 « « чзу(к)<t, х) , , ,
+тИ^<t,х) \ к = 1,2,...,Д.
2 г=1 ]■=1 йгг йг;.
3. Применение спектральной формы математического описания
Сведем задачу анализа к поиску коэффициентов разложения р . . функции ^, х) по ортонормированным функциям в
пространстве L2<T x Rn). Обозначим полную систему ортонор-мированных функций {е</0,/1,...,/п,t,х)}“. . =0, предполагая,
что они представляются в виде произведения
е</'о,*1,...,/„,t,х) = q<io,t) • р^,..., 1п,х), /'о,/1,...,/„ = 0,1,2,...,
причем функции {q<io,0}”=о и {р^,...,*п,х)}“..,п =о образуют
полные ортонормированные системы в пространствах L2<T) и L2<Rn) соответственно.
Применим спектральное преобразование [11, 12] к системе уравнений (2)-(5) с учетом условий (7), тогда
Р(п +1,п +1) -Ф(1)(п +1,0) - q(1,0;^) ® Ф0(п,0) =
= А(п +1, п +1) • Ф(1)(п +1,0) - (Л1 (п +1, п +1) +
+Л 2(п +1, п +1)) • Ф(1)(п +1,0) + Я1(п +1, п +1) •
• Ф( ^ )(п +1,0) + Н 2(п +1, п +1) • Ф( д) (п +1,0),
Р(п +1, п +1) • Ф(к) (п +1,0) = А(п +1, п +1) • Ф(к) (п +1,0) -(10) -Л1 (п +1, п +1) ^ )(п +1,0) + Л1(п +1, п +1) • Ф(к-1) (п +1,0), k = 2,..., Д1,
Р(п +1, п +1) • Ф( ^+1) (п +1,0) = А(п +1, п +1) • Ф( ^+1) (п +1,0) --Л 2 (п +1, п +1) • Ф( ^+1) (п +1,0) + Л 2(п +1, п +1) • Ф(1) (п +1,0), Р(п +1, п +1) • Ф№) (п +1,0) = А(п +1, п +1) • Ф№) (п +1,0) -
(12) -Л 2(п +1, п +1) -Ф^ )(п +1,0) + Л 2(п +1, п +1) • Ф№-1) (п +1,0), k = N. + 2,..., N.
В этих соотношениях Р(п + 1, п + 1) - спектральная характеристика оператора дифференцирования по времени с учетом значения функции в начальный момент; А(п + 1, п + 1), Ні(п + 1, п + 1) и Н2(п + 1, п + 1) - спектральные характеристики операторов А, Н\ и Н2, определенных выражениями (6) и (8); Лі(п + 1, п + 1) и Л2(п + 1, п + 1) - спектральные характеристики операторов умножения на функции А\(ґ) и Л2(0 соответственно; Фда(п + 1, 0) - спектральные характеристики функций ф(А%, х), k = 1,2,...,N. Все перечисленные спектральные характеристики определены относительно системы (е(/0,і1,...,іп,t,х)}“; г- =0. Далее, q(1, 0; t0) - матрица-столбец значений функций базисной системы ^(і0,0}“=0 в точке t0; Ф0(п, 0) - спектральная характеристика плотности вероятности ф0(х) начального состояния Х0, определенная относительно системы {р(^,...,іп,х)}“ г. =0.
Спектральная характеристика Ф(п + 1, 0) плотности вероятности ф(ї, х), называемая также обобщенной характеристической
функцией [11, 12], выражается следующим образом (Ф(п + 1, 0) - многомерная гиперстолбцовая матрица, образованная искомыми коэффициентами разложения (рі. і ):
N
(13) Ф(п +1,0) = ^Фда(п +1,0).
k=1
В основе соотношений (9)-(12) лежат определения, форма представления и свойства спектральных характеристик функций и линейных операторов, подробно изложенные в [12, 17, 18]. Определение спектральной характеристики, аналогичной Ні(п + 1, п + 1) и Н2(п + 1, п + 1), дано в [10, 12]. В этих же работах приведено представление спектральной характеристики А(п + 1, п + 1) с помощью спектральных характеристик операторов дифференцирования и умножения.
Система уравнений (9)-(12) - это система линейных матричных уравнений относительно неизвестных спектральных характеристик Ф®(п + 1, 0), или система линейных неоднородных алгебраических уравнений относительно элементов матриц Фда(п + 1, 0), т.е. коэффициентов разложения функций х) по функциям {е(і0,і1,...,іп,t,х)}“і . =0. Перейдем к решению
этой системы.
Из уравнений (10)-(12) следует, что
(Р(п +1, п +1) - А(п +1, п +1) + Л1(п +1, п +1)) • Ф(4г )(п +1,0) =
= Л1(п +1,п +1)• Ф(1с-1)(п +1,0), k = 2,...,Nl,
(Р(п +1, п +1) — А(п +1, п +1) + Л2 (п +1, п +1) ) •
• Ф( ^ -1) (п +1,0) = Л 2(п +1, п +1) • Ф(1) (п +1,0),
(Р(п +1, п +1) - А(п +1, п +1) + Л 2(п +1, п +1)) • Ф1^ )(п +1,0) =
= Л2(п +1,п +1)•Ф{ІІ-1)(п +1,0), k = N + 2,...,N,
т.е.
Ф№-1) (п +1,0) = Л-1 (п +1, п +1) •( Р(п +1, п +1) -
-А(п +1,п +1) + Л1(п +1,п +1))• Ф№)(п +1,0), k = 2,...,Nl,
Ф(1)(п +1,0) = Л 2'(п +1, п +1) • (Р(п +1, п +1) - А(п +1, п +1) + +Л2 (п +1, п +1) )• Ф( ^ -1) (п +1,0),
Ф(к-1) (п +1,0) = Л-1 (п +1, п +1) •( Р(п +1, п +1) -
- А(п +1, п +1) + Л 2(п +1, п +1) )• Ф(к} (п +1,0), k = N + 2,..., N, или кратко
Ф(к-1) (п +1,0) = Щ(п +1,п +1) • Ф(к) (п +1,0), k = 2,..., N1,
Ф(1) (п +1,0) = W2(n +1, п +1) • Ф( *+1 (п +1,0),
Ф(к-1)(п +1,0) = W2(n +1,п +1)-Ф^)(п +1,0), k = N1 + 2,...,N,
где
Щг(п +1, п +1) = Л-1 (п +1, п +1) • (Р(п +1, п +1) -—А(п +1, п +1) + Л^ (п +1, п +1)),
W2(n +1, п +1) = Л-1 (п +1, п +1) • (Р(п +1, п +1) -—А(п +1, п +1) + Л 2 (п +1, п +1)).
Таким образом,
Ф(1) (п +1,0) = -1 (п +1, п +1) • Ф(к) (п +1,0),
Ф(к)(п +1,0) = -к(п +1,п +1)^Ф(^(п +1,0), k = 1,...,NI,
Ф(1) (п +1,0) = - * (п +1, п +1) • Ф (к) (п +1,0),
Ф(к} (п +1,0) = Щ?-к (п +1, п +1) ^Ф(N](п +1,0), к = N +1,..., N, в частности
(14) Ф(1) (п +1,0) = Щ?1 -1 (п +1, п +1) •Ф(^(п +1,0),
(15) Ф(1)(п +1,0) = Щ^2 -1(п +1, п +1) ^Ф( N )(п +1,0).
Значит,
Щ*1 -1(п +1, п +1) • Ф( ^ )(п +1,0) =
= Щ^2-1 (п +1, п +1) • Ф( N) (п +1,0),
т.е.
Ф( ^(п +1,0) = Щ'1 -1 (п +1, п +1)]-1 •
• Щ^2-1 (п +1, п +1) • Ф( N) (п +1,0),
или
60
Ф( д) (п +1,0) = [ЖД2 1 (п +1, п +1)]-1 •
• -1 (п +1, п +1) •Ф(Ді)(п +1,0).
Перепишем уравнение (10) с учетом (14)-(16):
(л1 (п +1, п +1) • Жд (п +1, п +1) + Л 2(п +1, п +1) •
• ЖД1 1 (п +1, п +1) - Я1(п +1, п +1) - Н 2(п +1, п +1) •
• [ЖД2-1 (п +1, п +1)]-1 Ж*1 -1 (п +1, п +1)) • Ф( Д1)(п +1,0) =
= 9(1,0; ^) ®Ф0(п,0),
или
(л 2 (п +1, п +1) • Ж22 (п +1, п +1) + Л1(п +1, п +1) •
• Ж2Д 1 (п +1, п +1) - Н 2(п +1, п +1) - Н1(п +1, п +1) •
• [ЖД1 -1 (п +1, п +1)]-1 • Ж,*2 -1 (п +1, п +1) )• Ф( м) (п +1,0) =
= 9(1,0;^) ®Ф0(п,0), следовательно,
Ф(д?1)(п +1,0) = (л1(п +1, п +1) • Ж*1 (п +1, п +1) - Н1(п +1, п +1) + +{Л2 (п +1, п +1) - Н2(п +1, п +1) • [ЖД2 -1 (п +1, п +1)]-1} •
•ЖД1 -1(п +1, п +1))-1 •( q(1,0;to) ®Ф0(п,0)),
или
Ф(д )(п +1,0) = (л2 (п +1, п +1) •ЖД2(п +1, п +1) - Н2(п +1, п +1) + +{Л1 (п +1, п +1) - Н1(п +1, п +1) • [ЖД1 1 (п +1, п +1)]1} •
• ЖД2 -1 (п +1, п +1) )-1 •( 9(1,0; О ® Ф0 (п, 0)).
Выразим спектральную характеристику Ф(п + 1, 0), принимая во внимание (13):
N N
Ф(п +1,0) = ХФ(к)(п +1,0) + X ф(к)(п +1,0) =
к=1 к=N1 +1
N1
(19) = ^ЩЩ -к (п +1, п +1) •Ф^п +1,0) +
к=1 N2 -1
+ £ W2N2-к-1 (п +1, п +1) • Ф(N) (п +1,0).
к=1
В [10] было показано, что
X WN-к (п +1, п +1) = (Е (п +1, п +1) - WN (п +1, п +1) V
(20) к=1 ' ’
•(Е(п +1, п +1) - W(п +1, п +1))-', поэтому с учетом (17)
Ф(п +1,0) = [(Е (п +1, п +1) - Щ*1 (п +1, п +1) )•
• (Е (п +1, п +1) - Щ(п +1, п +1) )-1 +
+(Е (п +1, п +1) - Щ^2 -1(п +1, п +1) )•
• (Е (п +1, п +1) - Щ2(п +1, п +1) )-1 •
• [ ЩД2-1 (п +1, п +1)]-1 • Щ*1 -1 (п +1, п +1) ] • Ф( ^ )(п +1,0), или с учетом (18)
Ф(п +1,0) = [(Е (п +1, п +1) - ЩД1 (п +1, п +1)) •
• (Е (п +1, п +1) - Щ(п +1, п +1) )-1 • [Щ^1 -1(п +1, п +1)]-1 •
• Щ?2-1 (п +1, п +1) + (Е (п +1, п +1) - Щ?2-1 (п +1, п +1)) •
• (Е (п +1, п +1) - Ж,(п +1, п +1) )-1 Таким образом,
Ф( д) (п +1,0).
Ф(п +1,0) = [(Е (п +1, п +1) - ЩД1 (п +1, п +1)) •
•( Е (п +1, п +1) - Щ1(п +1, п +1)) 1 +(Е (п +1, п +1) --Щ^ -1 (п +1, п +1)) • (Е(п +1, п +1) - Щ2(п +1, п +1) )-1 •
(21) • [Щ^-1 (п +1, п +1)]-1 • Щ?1 -1 (п +1,п +1)] • (Л1 (п +1,п +1) •
• ЩД1(п +1, п +1) - Я1(п +1, п +1) + {Л 2(п +1, п +1) -
- Н2 (п +1, п +1) • [Щ^2-1 (п +1, п +1)]-1} •
• Щ^1 -1 (п +1, п +1) )-1 • (<?(1,0;^) ® Ф0(п,0)),
или
Ф(п +1,0) = [(Е (п +1, п +1) - 1 (п +1, п +1)) •
• (Е (п +1, п +1) - Щ(п +1, п +1) )-1 • Щ1 -1(п +1, п +1)]-1 •
• -1 (п +1, п +1) + (Е (п +1, п +1) - ЩД2 -1 (п +1, п +1)) •
(22) •(Е(п +1,п +1)-Щ2(п +1,п +1))-1 ] • (Л2(п +1,п +1) •
• Щ2щ(п +1, п +1) - Н 2(п +1, п +1) + {Л1(п +1, п +1) -- Н1 (п +1, п +1) • [Щ^1 -1 (п +1, п +1)]-1} •
• Щ^2-1 (п +1, п +1) )-1 • ( q(1,0; ^) ® Ф0 (п, 0)),
- решение рассматриваемой задачи в спектральной форме математического описания.
Для нахождения решения задачи анализа в пространстве функций времени и вектора состояния требуется применить формулу обращения [12]:
ад ад ад
, х) = X X... X Рад...п • ^ h,...,*п, t, х1 (t, х) е Т х ^.
г0 =0 . =0 .п =0
Обычно приближенно определяется конечное число коэффициентов р.. . , поскольку задача нахождения всех коэффициентов разложения не является тривиальной. В этом случае бесконечные матрицы в (9)-(13) заменяются конечными матрицами, тогда
и-1и -1 1„ -1
Ни х) * X X... X Юй...п • ^ il,..., и, t, хХ
г0 =0 . =0 .п =0
где натуральные числа Ь0, Ь1, ..., ип - выбранные порядки усечения спектральных характеристик, влияющие на точность решения. Проверка корректности расчетов обычно осуществляется путем сравнения результатов при использовании различных наборов базисных систем или усечений спектральных характеристик [18].
Для апробации разработанного метода рассмотрим стохастическую систему, заданную уравнением Ито
dX(t) = а sign X(t) dt + odW^) + dQ(t), X(0) = Х0, где X е R - состояние системы, t е [0,1]; а, а - числовые параметры; Щ(0 - скалярный стандартный винеровский процесс; Q(t) - скачкообразный процесс, определяемый интенсивностями А1, Х2, порядками N1, N2 и плотностями вероятности q1(t, у),
q2(t, у).
Для численных расчетов положим а = -0,5; а = 0,75; начальное состояние X0 имеет нормальное распределение с математическим ожиданием т0 = 1 и дисперсией D0 = 0,5; интенсивности и порядки заданы следующим образом: А1 = 6; Х2 = 1,5, N1 = N = 2; плотность распределения q1(t, у) приращения У1 - гауссовская с параметрами т1 = 0,5 и D1 = 0,1, q2(t, у) = = д(у + 0,2), т.е. величина приращения У2 постоянна и равна -0,2.
При решении задачи спектральным методом в качестве базисной системы ^(г0,t)}ад=0 были выбраны косинусоиды [12,
17, 18], а в качестве базисной системы {р(.1, х)}ад=0 - функции Эрмита [14] с параметрами т = 1 и D = 1, порядки усечения
и = и = 16.
Результаты вычислений для плотности вероятности состояния системы приведены на рис. 3, а на рис. 4 показаны ее сечения в моменты времени 0,0; 0,5; 1.0.
Рис. 3. Плотность вероятности состояния системы
Далее, на рис. 5 показано, как меняется плотность вероятности в конечном моменте времени в зависимости от порядка N1 (дополнительные расчеты проводились для случаев N1 = 1 и N = 3, при этом интенсивность Х1 определялась из условия X1/N1 = const, которое обеспечивает постоянство среднего числа скачков в единицу времени для эрланговского закона с параметрами Х1 и N1, т.е. N1 = 1, Х1 = 3 и N1 = 3, Х1 = 9). Видно, что максимальное значение для плотности вероятности и математическое ожидание состояния смещаются вправо. Это объясняется ростом вероятности выбора приращения Y1, а не Y2: при Х1 = 3 эта вероятность равна 3/(3 + 1,5) = 2 / 3, при Х1 = 6 -
6/(6 + 1,5) = 4/5 и при Х1 = 9 - 9/(9 + 1,5) = 6/7.
Асимметрию можно объяснить положительной средней величиной приращения Y1, которая по модулю превосходит детерминированное приращение Y2, и, конечно, соотношением интенсивностей Х1, Х2.
Наконец, на рис. 6 показано, как зависит плотность вероятности в конечном моменте времени от различных усечений спектральных характеристик (дополнительные расчеты проводились для случаев L0 = L1 = 8 и L0 = L1 = 24 при номинальных параметрах N1 = 2, Х1 = 6).
Как следует из графиков, при выбранных порядках усечения (8; 16; 24) разница в приближенных решениях не столь заметна, однако при другом масштабе видно (см. рис. 7), что с ростом усечения растет и точность. В первую очередь это проявляется в минимизации интервалов, где оценка плотности вероятности отрицательна (появление таких интервалов при приближенном решении неизбежно и вызвано свойствами выбранных базисных функций).
Г = 0,0 ^ Г = 0,5 X*
Лч / Г = 1,0
4 -2 0 2 4 6
Рис. 4. Сечения плотности вероятности состояния системы в различные моменты времени
^=1- =3
/
/ N^=2
4 -2 0 2 4 6
Рис. 5. Сечения плотности вероятности состояния системы в конечном моменте времени при различных значениях параметров Х1 и N1
Рис. 6. Сечения плотности вероятности состояния системы в конечном моменте времени при различных порядках
усечения
Отметим, что все численные расчеты были выполнены с помощью специализированного программного обеспечения Spectrum [15].
0.2
0.1
/
16 24 =U-« /
-3 -2 -1 0
Рис. 7. Сечения плотности вероятности состояния системы в конечном моменте времени при различных порядках усечения
4. Заключение
В работе рассмотрено применение спектральной формы математического описания к задаче вероятностного анализа стохастических систем, которые характеризуются наличием разрывов (скачков) траекторий, образующих гиперэрланговский поток событий. Получены соотношения для нахождения плотности вероятности вектора состояния в спектральной форме математического описания систем управления. Использование гиперэр-ланговских потоков дает возможность учитывать более сложный характер поведения траекторий вектора состояния, а применение спектральной формы математического описания позволяет свести систему интегро-дифференциальных уравнений к системе линейных алгебраических уравнений относительно коэффициентов разложения искомой плотности вероятности по некоторой ортонормированной системе функций.
Рассмотренный подход упрощает процесс решения задачи анализа, делая его удобным для применения современных высокопроизводительных вычислительных систем.
Литература
1. АВЕРИНА Т.А. Модифицированный алгоритм статистического моделирования систем со случайным периодом квантования // Вестник СГТУ. - 2011. - №4(62), Вып. 4. -С.212-218.
2. АВЕРИНА Т А., РЫБАКОВ К.А. Новые методы анализа
воздействия пуассоновских дельта-импульсов в задачах радиотехники // Журнал радиоэлектроники. - 2013. - №1. -[Электронный ресурс]. - иЯЬ: http://jre.cplire.ru/jre/
contents.html (дата обращения: 09.09.2013).
3. АНУЛОВА СВ., ВЕРЕТЕННИКОВ А.Ю., КРЫЛОВ Н.В. И ДР. Стохастическое исчисление // Итоги науки и техники. Современные проблемы математики. Фундаментальные направления. Т. 45. - М.: ВИНИТИ, 1989. - С. 42-79.
4. АРТЕМЬЕВ В.М., ИВАНОВСКИЙ А.В. Дискретные системы управления со случайным периодом квантования. -М.: Энергоатомиздат, 1986. - 96 с.
5. ВЕРЕТЕННИКОВ А.Ю. О сильных решениях стохастических уравнений Ито со скачками // Теория вероятностей и ее применения. - 1987. - Т. 32., Вып. 1. - С. 159-163.
6. ГОРИЦКИЙ Ю .А., КАЗАКОВ В.А. Дискретизация случайных процессов с конечным множеством состояний и эрлан-говским временем пребывания // Известия РАН. Теория и системы управления. - 2011. - №6. - С. 14-27.
7. КАЗАКОВ И.Е., АРТЕМЬЕВ В.М., БУХАЛЕВ В.А. Анализ систем случайной структуры. - М.: Физматлит, 1993. -272 с.
8. КОЖЕВНИКОВ А.С. Программное обеспечение для стати-
стического моделирования и анализа случайных процессов со скачками, описывающих динамику цен акций предприятий авиационной отрасли // Труды МАИ. - 2012. - № 59.
- [Электронный ресурс]. - иЯЬ:
http://www.mai.ru/science/trudy (дата обращения: 09.09.2013).
9. КОЖЕВНИКОВ А.С. Математические модели динамики цены акций с гиперэрланговскими скачками // Научный альманах. Вып. 17: Материалы IX научно-практической конференции молодых ученых и студентов «Инновационный менеджмент в аэрокосмической промышленности». - М.: Изд-во «Доброе слово», 2013. - С. 180-186.
10. КОЖЕВНИКОВ А.С., РЫБАКОВ К.А. Анализ нелинейных стохастических систем управления с импульсными воздействиями, образующими эрланговские потоки событий // Научный вестник МГТУ ГА. - 2012. -№184(10). - С. 37-45.
11. ПАНТЕЛЕЕВ А.В., РЫБАКОВ К.А. Анализ нелинейных стохастических систем управления в классе обобщенных характеристических функций // Автоматика и телемеханика. - 2011. - №2. - С. 183-194.
12. ПАНТЕЛЕЕВ А.В., РЫБАКОВ К.А., СОТСКОВА И.Л.
Спектральный метод анализа нелинейных стохастических систем управления. - М.: Вузовская книга, 2006. - 392 с.
13. ПУГАЧЕВ В.С., СИНИЦЫН И.Н. Стохастические дифференциальные системы. Анализ и фильтрация. - М.: Наука, І990. - 632 с.
14. РОМАНОВ В.А., РЫБАКОВ К.А. Спектральные характе-
ристики операторов умножения, дифференцирования и интегрирования в базисе обобщенных функций Эрмита // Труды МАИ. - 20і0. - №39. - [Электронный ресурс]. -URL: http://www.mai.ru/science/trudy (дата обращения:
09.09.2013)
15. РЫБАКОВ К.А. Программное обеспечение спектрального
метода Spectrum // Труды МАИ. - 2003. - №І4. - [Электронный ресурс]. - URL: http://www.mai.ru/science/
trudy (дата обращения: 09.09.20і3).
16. РЫБАКОВ К.А. Спектральные характеристики операторов умножения, дифференцирования и интегрирования в базисе обобщенных функций Лагерра // Дифференциальные уравнения и процессы управления. - 20І2. - №і. - С. ІІ4-і4і. - [Электронный ресурс]. - URL: http://www.math.spbu.ru/diffjournal (дата обращения:
09.09.2013).
17. РЫБИН В.В. Моделирование нестационарных непрерывнодискретных систем управления спектральным методом в системах компьютерной математики. - М.: Изд-во МАИ, 20ІІ. - 220 с.
18. СОЛОДОВНИКОВ В.В., СЕМЕНОВ В.В. Спектральная теория нестационарных систем управления. - М.: Наука, І974. - 336 с.
19. ТИХОНОВ В.И., МИРОНОВ М.А. Марковские процессы. -М.: Советское радио, І977. - 488 с.
20. ФЕДОСОВ Е.А., ИНСАРОВ ВВ., СЕЛИВОХИН ОС.
Системы управления конечным положением в условиях противодействия среды. - М.: Наука, І989. - 272 с.
21. CONT R., TANKOV P. Financial Modelling with Jump Processes. - Chapman and Hall, 2004. - 560 p.
22. IWANKIEWICZ R., NIELSEN S.R.K. Advanced Methods in Stochastic Dynamics of Non-Linear Systems. - Aalborg tekniske Universitetsforlag, 1999. - 276 p.
23. HANSON F.B. Applied Stochastic Processes and Control for Jump-Diffusions: Modeling, Analysis, and Computation. -SIAM, 2007. - 472 p.
24. 0KSENDAL B., SULEM A. Applied Stochastic Control of Jump Diffusions. - Springer, 2007. - 266 p.
25. TELLIER M., IWANKIEWICZ R. Response of linear dynamic systems to non-Erlang renewal impulses: stochastic equation approach // Probabilistic Engineering Mechanics. - 2005. -Vol. 20, №4. - P. 281-295.
SPECTRAL METHOD FOR STOCHASTIC SYSTEMS WITH DISCONTINUOUS TRAJECTORIES DESCRIBED BY RANDOM MIXTURE OF ERLANG DISTRIBUTIONS Alexander Kozhevnikov, Moscow aviation institute, graduate student (exequit@yandex.ru).
Konstantin Rybakov, Moscow aviation institute, Cand.Sc., assistant professor (rkoffice@mail.ru).
Abstract: We consider stochastic control systems with impulses generated by hyper-Erlang flows of events leading to jumps in system trajectories. We employ the spectral form of mathematical description of the system to find the probability density function of its states.
Keywords: analysis problem, hyper-Erlang distribution, jump process, spectral method, stochastic system.
Статья представлена к публикации членом редакционной коллегии А.П. Курдюковым
Поступила в редакцию 17.01.2013.
Опубликована 30.09.2013.
