CC BY
61
УДК 519.676
ПОСТРОЕНИЕ И ОБОСНОВАНИЕ СТАТИСТИЧЕСКИХ АЛГОРИТМОВ МОДЕЛИРОВАНИЯ РЕШЕНИЯ СИСТЕМ СО СЛУЧАЙНОЙ СТРУКТУРОЙ, ЗАДАННОЙ СТОХАСТИЧЕСКИМИ ДИФФЕРЕНЦИАЛЬНЫМИ УРАВНЕНИЯМИ
© Т.А. Аверина
Ключевые слова: статистическое моделирование; системы со случайной структурой; стохастические дифференциальные уравнения; пуассоновский поток; метод «максимального сечения».
В данной работе построены алгоритмы статистического моделирования решения систем со случайной структурой с независимыми распределенными переходами. Исследована слабая сходимость численных методов, построенных на основе разработанных алгоритмов.
Под системами со случайной структурой понимаются динамические системы, поведение которых на случайных интервалах времени характеризуется различными структурами и описывается различными уравнениями [1].
Состояние системы со случайной структурой характеризуется смешанным процессом [y(t),s(t)]T , где s(t) - дискретный случайный процесс с конечным множеством состояний {1, 2,...,S} , S - число структур системы, а y(t) - d -мерный непрерывный случайный процесс, описываемый при условии s(t) = l стохастичегким дифференциальным уравнением (СДУ) в форме Стратоновича:
dy(t) = a(i) (t, y(t))dt + a(i) (t, y(t))dw(t), y(to) = yo. (1)
Здесь t € [to,T] ; w(t) - m -мерный стандартный винеровский процесс, не зависящий от yo ; a(i)(t, y) : [to, T] x Rd ^ Rd - вектор-функция размера d ; a(i)(t, y) : [to,T] x Rd ^ Rdxm -матричная функция размера d x m ; l - номер структуры, l = 1, 2,..., S .
Вероятность перехода дискретного случайного процесса s(t) удовлетворяет условию [1]:
P(s(t + At) = r | s(t) = l, y(t) = y) = Vir(t)At + o(At),
P(s(t + At) = l | s(t) = l, y(t) = y) =1 - vi(t)At + o(At),
s(to) = so, l,r = 1,2,..., S, l = r, (2)
где функция Vir(t) : [to,T] x Rn ^ [0, называется интенсивностью перехода из l-й
структуры к r -й, причем Vi(t) = S=i=i vir(t) . Известно, что ненормированные плотности распределения для каждой структуры удовлетворяют системе обобщенных уравнений Фоккера-Планка-Колмогорова [1].
Условие (2) обеспечивает отсутствие нескольких переключений процесса s(t) за малый интервал времени At. Такое поглощение реализаций называют пуассоновским [1]. Моменты смены структуры представляют пуассоновский поток. Поэтому алгоритмы статистического моделирования пуассоновских точечных полей, предложенные в работах [2, 3], можно использовать для моделирования моментов смены структуры. Статистические алгоритмы моделирования решения систем со случайной структурой c постоянными интенсивностя-ми перехода можно найти в работах [4, 5]. В данной работе мы построим статистические алгоритмы для переменных интенсивностей перехода.
Напомним, что обобщенное экспоненциальное распределение случайной величины £ с функцией интенсивности Л имеет вид:
p(x) = Л(х)е- J0 x(u)du, x > 0, Л(х) ^ 0 Л(и)^и =
Jo
Тогда при условии s(t) = l промежутки rir времен перехода из l -й структуры к r -й имеют обобщенное экспоненциальное распределение с плотностями:
pTir (x) = Vlr (t + x)e^ox (i+u)du, x> 0, l, r = 1, 2, r = l, (3)
и смена структуры происходит в момент времени t = t + то , где
то = min Tir. (4)
r=l
Теорема 1. Пусть задана система со случайной структурой (1) с независимыми распределенными переходами (2) и Tir (r = 1,2,...,S, r = l) - промежутки времен перехода из l -й к r -й структуре, распределенные при условии s(t) = l с плотностями (3). Тогда, при условии s(t) = l :
1) случайная величина т0 (4) имеет обобщенное экспоненциальное распределение c плотностью
Рт0(т) = Vi(t + т)e^oT Vi(t+u)du, т > 0; (5)
2) номер новой структуры s(t + т0) является случайной величиной п, для которой
P(п = г|то = т) = Vlr(t + т), r = 1,2,..., S, r = l. (6)
Vl(t + т)
Доказательство теоремы основано на свойствах экспоненциального распределения. В случае постоянных интенсивностей перехода
S
Vlr (t) = Vir, Vi = Vi = ^ Vir
r=1=l
моделирование времени нахождения системы в текущей структуре осуществляется по формуле
т = —ln(a)/Vl.
В общем случае для моделирования распределения (5) можно использовать метод « максимального сечения > с постоянной
S
Vlr (t) ^ Vir, Vi (t) ^ Vi = ^ Vlr, t > 0
r=1=l
или переменной мажорантой
S
Vlr (t) < Vlr (t), Vi(t) < Vi(t) = Vlr (t), t > 0, (7)
i=i
если функции Vi(t) вычисляются достаточно просто.
Обозначим tk - момент смены структуры и Sk - номер новой структуры, тогда
s(t) = sk при tk ^ t < tk+1, tK = T, k = 0,1,...,K.
Используя способы моделирования обобщенного экспоненциального распределения с плотностью (5), рассмотренные в работе [6], можно сформулировать следующие алгоритмы моделирования процесса смены структуры для решения систем со случайной структурой (1) с независимыми распределенными переходами (2).
Алгоритм 1а ( для систем с постоянными интенсивностями перехода):
0) k:=0; моделируем Sk согласно заданному начальному условию S(t0) ;
1) l := Sk ; а := rand; т = — ^ (время нахождения в l -й структуре (4));
2) tk+1 = tk + т ; если tk+1 ^ T, то STOP;
3) а := rand;
r-1_ r _
— < а ^ У, то Sk+1 = r;
i=1=l Vl i=1=l Vl
4) k := k + 1, идем на 1).
Алгоритм 1b ( использование метода «максимального сечения»):
0) k:=0; моделируем Sk согласно заданному начальному условию s(t0);
1) l := Sk ; моделируем т - возможное время нахождения системы в текущей структуре по методу «максимального сечения»:
1.1) т := 0 ;
1.2) моделируем £ как значение случайной величины с условной плотностью p(x) = = V (tk + x)e- ^(tfc +u)du, x > 0 и полагаем т := т + £ ;
1.3) tk+1 = tk + т ; если tk+1 ^ T, то STOP;
1.4) а := rand; если а > v'(tfc+l) , то идем на 1.2);
' Vl(tk +1) 7 >'
2) а := rand ;
Evli(tk+1) „ . V^ vli(tk+1)
----<а ^ > -—--, то Sk+1 = r;
i=1=l vl(tk+1) i=1=l vl(tk+1)
3) k := k + 1, идем на 1).
Метод «максимального сечения» основан на правиле прореживания пуассоновского потока [7]: оставляя точки пуассоновского потока интенсивности V(t) с вероятностью
= Щ • (8)
мы получаем пуассоновский поток интенсивности
Vl(t) ■ p(t) = V(t) ■ = Vl(t).
В некоторых случаях, при большом числе структур S, достаточно трудоемко вычислять сумму функций при вычислении вероятностей прореживания (8). В этих случаях для уменьшения трудоемкости можно использовать представление вероятности в виде
p(t) = f V# = f ■ VlM. (9)
Vl(t) i=lVl(t) (t) Vli(t)
Метод «максимального сечения», использующий представление (9), называют методом «мажорантной частоты» [8].
Алгоритм 1с ( использование метода «мажорантной частоты»):
0) k:=0; моделируем Sk согласно заданному начальному условию s(t0);
1) l := Sk; моделируем возможное время перехода из l -й структуры и возможный номер новой структуры по методу «мажорантной частоты»:
1.1) т := 0 ;
1.2) моделируем £ как значение случайной величины с условной плотностью p(x) = = Vl(tk + x)e- Jo vi(tk +u)du, x > 0 и полагаем т := т + £ ;
1.3) полагаем tk+i = tk + т ; если tk+i ^ T, то STOP;
1.4) а := rand;
Ej Vii(tk+i) „ . j Vii(tk+i)
—-- < а ^ > —--, то r = j;
i=i=i Vi(tk+i) Vi(tk+i) ' ''
1.5) а := rand; если а > v}r(tk+1\ , то идем 1.2);
' ' vlr(tk +1) 7 "
2) Sk+i = r; k := k + 1, идем на 1).
Теорема 2. 1) Целочисленный случайный процесс s(t) , построенный с помощью алгоритма 1а, является процессом смены структуры для системы (1), (2) с постоянными интенсивностями перехода {Vir = Vir} , а целочисленные случайные процессы, построенные с помощью алгоритмов 1b, 1c являются процессами смены структуры для системы (1), (2) с переменными интенсивностями перехода {Vir(t)} , при условии (7) .
2) Построенные алгоритмы реализуется с конечной трудоемкостью.
Доказательство теоремы 2 основано на теоремах, доказанных в работах [3, 8].
Замечание. В работе [9] был предложен экономичный способ моделирования случайных величин, плотности распределения которых представляют собой взвешенные суммы (смеси) эффективно моделируемых плотностей. В работе [3] эта методика формулируется и обосновывается в многоэтапном варианте и применяется для моделирования неоднородных пуассоновских ансамблей с использованием последовательности «исключений» по одному случайному числу. На основании теоремы, доказанной в [3], можно модифицировать алгоритмы 1b и 1c с использованием моделирования только одного случайного числа а.
Мы построили статистические алгоритмы моделирования процесса смены структуры, то есть совокупности пар точек {tk, Sk} , задающих момент смены структуры tk и номер новой структуры Sk . Численные методы решения задачи Коши для стохастических дифференциальных уравнений (1) можно найти в [10,11] . С учетом построенных алгоритмов и численных методов решения СДУ, можно записать следующий
Алгоритм2 (статистического моделирования решения системы со случайной структурой с независимыми распределенными переходами (1) , (2) ):
1) моделируем процесс смены структуры {tk, Sk} , k = 0, ...K , sk = T одним из алгоритмов, построенных в предыдущем разделе;
2) k:=0;
3) на интервале [tk, tk+i] численным методом решения СДУ решаем уравнения (1) для Sk -й структуры, находим yk+i - вектор состояния системы в момент времени tk+i ;
4) если k + 1 < K , то
4а) довычисляем yk+i согласно qsksk+1 (y,t|y',t) - условной функции плотности вероятности восстановления Sk+i -й структуры из Sk -й;
4b) k := k + 1 и идем на 3);
5) если k + 1 = K , то STOP.
Вид функции qir определяется физическим содержанием задачи и характеризует начальные условия при восстановлении процесса в r -й структуре при переходе из l -й структуры. Могут иметь место различные условия восстановления, определяемые функциями qir . Рассмотрим случай точного восстановления
qir (y,t|y',t) = 5(y - y'),
т. е. конечные условия процесса в l -й структуре совпадают с начальными в r -й структуре.
Выбор численного метода решения конкретной системы СДУ и шага интегрирования h определяются видом этой системы и требуемой точностью вычисления вероятностных характеристик выходных процессов. Для различных структур могут использоваться разные численные методы с различными шагами интегрирования.
Пусть для решения l -й структуры используется численный метод, который слабо сходится с порядком pi на решении задачи Коши для СДУ, и временная сетка {Tn} , n = = 0,1, ...N включает равномерную сетку с шагом h и моменты скачков {tk} .
Определение. Численный метод решения систем со случайной структурой слабо сходится с порядком p, если для любой достаточно гладкой функции f (y)
max
1<ra<N
E(f (y(in))/y(io) = yo) - E(f (yfc)/y(to) = yo)| = O(hp), h ^ 0
где yn - численное решение в момент времени Tn , hn = Tn+i — Tn , h = max hn - мак-
l^ra^N
симальный шаг интегрирования, Tn = T, E - операция вычисления математического ожидания.
Теорема 3. Численный метод решения системы со случайной структурой (1) с распределенными переходами (2), заданный алгоритмом 2, имеет р-й порядок слабой сходимости:
p = min pi,
где pi - порядок слабой сходимости численного метода решения задачи Коши для СДУ, используемый для решения l -й структуры, l = 1, ...S .
Доказательство теоремы 3 основано на использовании формулы полной вероятности и теоремы сходимости, доказанной в [10].
ЛИТЕРАТУРА
1. Казаков И.Е., Артемьев В.М., Бухалев В.А. Анализ систем случайной структуры. М.: Физматлит, 1993.
2. Михайлов Г.А., Аверина Т.А. Алгоритм < максимального сечения » в методе Монте-Карло // ДАН. 2009. Т. 428. № 2. С. 163-165.
3. Аверина Т.А. Новые алгоритмы статистического моделирования неоднородных пуассоновских ансамблей // Журн. вычисл. матем. и матем. физ. 2010. Т. 50. № 1. С. 16-23.
4. Аверина Т.А. Численный анализ систем управления динамическими объектами со случайными изменениями структуры // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2011. Т. 16. Вып. 4. С. 1023-1026.
5. Аверина Т.А. Исследование влияния пуассоновских дельта-импульсов в задачах радиотехники // Вестник Тамбовского университета. Серия Естественные и технические науки. Тамбов, 2013. Т. 18. Вып. 5. С. 2431-2433.
6. Михайлов Г.А., Войтишек А.В. Численное статистическое моделирование. Методы Монте-Карло. М.: Изд. центр < Академия » , 2006.
7. Беляев Ю.К. Пуассоновский точечный процесс // Вероятность и математической статистика. Энциклопедия. М.: Научное изд-во, 1999. С. 525-526.
8. Михайлов Г.А., Рогазинский С.В. Модифицированный метод < мажорантной частоты » для численного моделирования обобщенного экспоненциального распределения // ДАН. 2012. Т. 444. № 1. С. 28-30.
9. Михайлов Г.А. К вопросу о построении экономичных алгоритмов моделированимя случайных величин // Журн. вычисл. матем. и матем. физ. 1966. Т. 6. № 6. С. 1134-1136.
10. Artemiev S.S., Averina T.A. Numerical Analysis of Systems of Ordinary and Stochastic Differential Equations. Utrecht, the Netherlands: VSP, 1997.
11. Аверина Т.А. Устойчивые численные методы решения стохастических дифференциальных уравнений в смысле Стратоновича // Вестник Бурятского Гос. Унив. 2012. № 9. С. 91-94.
БЛАГОДАРНОСТИ: Работа поддержана грантом РФФИ (проект № 14-01-00787) и грантом «Научные школы» НШ-5111.2014.1.
Поступила в редакцию 11 июня 2015 г.
Averina T.A. CONSTRUCTION AND JUSTIFICATION OF STATISTICAL ALGORITHMS FOR SIMULATION OF SWITCHING DIFFUSION, GIVEN BY STOCHASTIC DIFFERENTIAL EQUATIONS
The statistical algorithms for simulation of switching diffusion with independent distributed transitions are constructed. Weak convergence of these methods is investigated.
Key words: statistical simulation; switching diffusion; stochastic differential equations; Poisson flow; maximum cross section method.
Аверина Татьяна Александровна, Институт вычислительной математики и математической геофизики СО РАН, Новосибирский государственный университет, г. Новосибирск, Российская Федерация, кандидат физико-математических наук, старший научный сотрудник ИВМиМГ СО РАН, доцент кафедры вычислительной математики НГУ, e-mail: ata@osmf.sscc.ru
Averina Tatyana Aleksandrovna, Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk State University, Novosibirsk, the Russian Federation, Candidate of Physics and Mathematics, Senior Researcher of ICM&MG SB RAS, Associate Professor of the Computational Mathematics Department of NSU, e-mail: ata@osmf.sscc.ru
