ВЫЧИСЛЕНИЕ СРЕДНЕЙ СКОРОСТИ РОСТА ВЕКТОРА СОСТОЯНИЙ СТОХАСТИЧЕСКОЙ СИСТЕМЫ С СИНХРОНИЗАЦИЕЙ СОБЫТИЙ*
Н. К. Кривулин
С.-Петербургский государственный университет, д-р физ.-мат. наук, доцент, [email protected]
1. Введение. Динамические модели, которые описывают реальные системы в технике, экономике, управлении и других областях, нередко включают различные механизмы синхронизации событий, происходящих в системе. Например, в производственной системе процесс изготовления новой партии продукции на предприятии начинается не раньше, чем завершится поступление на это предприятие новой партии сырья. В системе передачи данных один узел начинает посылать очередное сообщение другому узлу, как только будет получено подтверждение об успешном приеме предыдущего сообщения, а в системе электронной торговли транспортировка товара заказчику начинается, как только будет осуществлено списание денежных средств со счета заказчика.
Во многих случаях динамика таких систем может быть описана при помощи векторных уравнений вида
г(к) = Л(к) <8> г (к — 1),
которые являются линейными в смысле некоторого идемпотентного полукольца [1—4]. Особый интерес представляют стохастические модели систем с синхронизацией, в которых матрица уравнения Л(к) является случайной. При этом целью исследования часто является определение средней асимптотической скорости роста вектора состояний системы г (к) при к ^ ж, которая обычно имеет смысл среднего времени одного цикла работы системы (производственного цикла, сеанса связи, банковской транзакции и т. п.). Величина, обратная среднему времени цикла, может рассматриваться как пропускная способность системы.
Еще одной актуальной задачей, тесно связанной с вычислением средней скорости роста вектора состояний, является нахождение предельных вероятностей, которые описывают поведение системы при к ^ ж. Например, при анализе некоторых систем может потребоваться вычисление вероятности совпадения моментов времени наступления различных событий или вероятности того, что интервал времени между событиями равен некоторому фиксированному числу.
Вычисление средней скорости роста вектора состояний может оказаться довольно трудной задачей даже для простых моделей систем. Большинство известных результатов получены для систем с матрицей второго порядка и экспоненциальным распределением ее элементов [4-8].
Рассмотрим матрицу
* Работа выполнена при финансовой поддержке РФФИ (грант №09-01-00808). © Н. К. Кривулин, 2011
Пусть элементы матрицы образуют последовательности {ак}, {вк}, {7к} и }, которые состоят из независимых случайных величин и сами являются независимыми. Элементы каждой последовательности имеют общее экспоненциальное распределение с параметрами V, а и т. Вычисление средней скорости роста вектора состояний для такой системы может быть осуществлено на основе алгебраического подхода в [4], при котором задача сводится к нахождению решения линейной системы уравнений и вычислению некоторого линейного функционала от полученного решения.
Если часть элементов матрицы системы является неотрицательными константами, то решение может быть найдено в явном виде. В частности, для систем с матрицей с нулевыми константами в работе [4] получены следующие результаты (далее указывается вид матрицы А(к) и выражение для средней скорости роста Л):
Известные решения для систем с матрицей, среди элементов которой есть произвольная константа с > 0, включают следующие случаи [7, 8]:
В настоящей работе рассматривается задача вычисления средней скорости роста вектора состояний для системы с матрицей второго порядка, у которой один элемент, стоящий на диагонали, является случайной величиной с экспоненциальным распределением, а все другие элементы равны одной и той же неотрицательной константе.
Решение задачи включает замену переменных, в результате которой вместо случайных координат вектора состояний вводятся новые случайные величины, анализ которых оказывается более удобным. Затем осуществляется построение и исследование сходимости соответствующих этим величинам последовательностей одномерных функций распределения. Средняя скорость роста вычисляется как среднее значение предельного распределения. В качестве дополнительного результата, который представляет самостоятельный интерес, получены выражения для предельных вероятностей некоторых событий в системе, включая вероятность равенства координат вектора состояний.
2. Стохастическая динамическая система. Рассмотрим стохастическую динамическую систему, эволюция которой для всех к — 1, 2,... определяется векторным уравнением
где z(k) —вектор состояний, А(к) —переходная матрица системы. Умножение матрицы на вектор понимается в смысле идемпотентного полукольца со скалярными операциями вычисления максимума в роли сложения ® и арифметического сложения в роли умножения 0.
^т (^ + т )(р? + т2 ) 4v2 + 7 vа + 4а2
Л — с +
Л — с +
z(k) — А(к) 0 z(k — 1),
Зададим матрицу и вектор состояний системы в виде
А(к)=(а с )■ =(х;к>). z(o)=(о
где с — неотрицательная константа, {«к} — последовательность независимых случайных величин, которые имеют общее экспоненциальное распределение с функциями распределения и плотности
„ , . (0, если Ь < 0, „ , , (0, если Ь < 0,
Ра(^) = N /а(Ь) = N
11 — е если Ь> 0; 1^е , если Ь > 0.
Векторное уравнение можно представить с помощью скалярных уравнений в идем-потентном полукольце
ж(к) = ак 0 ж(к — 1) © с 0 у (к — 1), у (к) = с 0 ж(к — 1) © с 0 у (к — 1),
или с помощью обычных уравнений
ж(к) = шах(ж(к — 1) + ак, у(к — 1) + с), у(к) = шах(ж(к — 1) + с, у(к — 1) + с).
3. Средняя скорость роста вектора состояний. Средняя (асимптотическая) скорость роста вектора состояний системы определяется в терминах полукольца как предел
Л = Иш |Ик)||1/к, к—
который с использованием обычных обозначений записывается в виде
А= Ит — тах(ж (к), у (к)), к—к
С помощью эргодической теоремы из [4] нетрудно показать, что для рассматриваемой системы указанный предел существует с вероятностью 1, а также выполняется равенство
Ит — Е тах(ж(/г), у (к)) = А. к—к
В некоторых случаях среднюю скорость роста вектора состояний системы можно вычислить без особого труда. Предположим, что с = 0. Тогда выполняется неравенство
ж(к) = шах(ж(к — 1) + ак, у(к — 1)) > шах(ж(к — 1), у(к — 1)) = у(к), а вычисление средней скорости роста сводится к нахождению величины
Л= Ит —Етах(ж (к),у(к)) = Ит —Еж (к), к—к к—к
Кроме того, из первого уравнения получаем
ж(к) = ж(к — 1) + «к = • • • = «1 +--+ «к,
откуда прямо следует, что
Л = Е«1 = 1/^.
Решение задачи для с > 0 требует дополнительных вычислений, которые будут представлены в следующих разделах.
4. Преобразование и предварительный анализ модели. Рассмотрим общий случай, когда с — произвольная неотрицательная константа. Сначала выполним замену переменных, предложенную в работе [5], которая позволяет свести решение задачи к анализу сходимости последовательностей одномерных функций распределения вероятностей.
Введем новые переменные состояния:
X (к) = ж(к) — ж(к — 1),
У (к) = у (к) — ж(к).
Имеем систему уравнений
X (к) = шах(ак, У (к — 1) + с),
У (к) = шах(0, У (к — 1)) — шах(ак — с, У (к — 1)).
Легко видеть, что выполняются неравенства
X (к) > «к > 0, У (к) < с — «к < с.
Кроме того, ясно, что
ж(к) = X (1) + ••• + X (к).
Определим функции распределения
Фк(Ь) = P{X(к) < Ь}, Фк(Ь) = Р{У(к) < Ь}.
Из первого уравнения системы следует равенство
Фк(Ь) = Р{шах(«к,У(к — 1) + с) < Ь} = ^а(Ь)Фк_ 1 (Ь — с).
Нетрудно проверить, что имеется положительная вероятность
Цк = Р{У(к) =0} = Р{У(к — 1) > шах(0, «к — с)} > 0.
Тогда с учетом предыдущего равенства для всех к > 1 выполняется
Рк = P{X(к) = с} > 0.
Заметим, что рк обозначает вероятность того, что приращение первой координаты вектора состояний при завершении цикла к будет равно с, а величина Цк — вероятность того, что обе координаты вектора z(k) будут совпадать.
Ясно, что константа с > 0 для X (к) и ноль для У (к) являются единственными значениями, которые эти случайные величины принимают с положительной вероятностью.
5. Сходимость последовательностей распределений. Исследуем последовательность функций распределения {Фк}. Заметим, что
. . (О, если Ь < 0, Г1 — ^а(с — Ь), если Ь < 0,
Ф°(Ь) = -п ф1(Ь)^ а “п
I 1, если Ь > 0, 11, если Ь > 0.
По формуле полной вероятности для всех к > 1 запишем
/* оо
Фк (Ь) = / Р{У (к) <Ь|«к = и}/а(и)^и.
°
Найдем условную вероятность
Р{У (к) < Ь|«к = и} = Р{шах(0, У (к — 1)) — шах(и — с, У (к — 1)) < Ь} =
'0, если и < с — Ь, Ь < 0,
Фк_ 1 (и — с + Ь), если и > с — Ь, Ь < 0,
1 — Фк-1 (—Ь), если и < с — Ь, Ь > 0,
_ 1, если и > с — Ь, Ь > 0.
Подстановка в формулу полной вероятности дает рекуррентное уравнение
ОО
Фк_ 1(и)/а(и + с — Ь)йи, если Ь < 0,
Фк(Ь) ^ 1 — Фк-1(—Ь)^а(с — Ь), если 0 < Ь < с,
если Ь > с.
Учитывая, что при 0 < Ь < с выполняется равенство
/* О
Фк (—Ь) = / Фк_1 (и)/а(и + с + Ь)йи,
°
рекуррентное уравнение можно записать в виде
• оо
Фк_ 1(и)/а(и + с — Ь)йи, если Ь < 0,
Фк(Ь) = <
а V
°
оо
1 — ^а (с — Ь) / Фк_2 (и)/а(и + с + Ь)йи, если 0 < Ь < с,
°
1, если Ь > с.
Подставляя выражение для экспоненциального распределения, получаем
Фк (Ь) = <
^е м(с ^ Фк_1 (и)е м“йи, если Ь < 0,
°
ОО
1 — ^е_м(с+4)(1 — е_м(с_4)) / Фк_2(и)е_м“йи, если 0 < Ь < с,
°
1, если Ь > с.
Определим числовую последовательность {ак} с элементами
/* О
а° = а1 = 1, ак = ^ Фк(и)е_м“йи, к > 1.
°
С учетом обозначения
С = е-^с
получим
{Сак_1ем*, если Ь < 0,
1 — Сак_2 (1 — Сем*)е_м*, если 0 < Ь < с,
1, если Ь > с.
Очевидно, что имеется взаимно однозначное соответствие между последовательностью функций {Фк} и числовой последовательностью {ак}.
Составим рекуррентное уравнение для числовой последовательности
ак = 1 — 2^(1 — С')‘2а1~-2-
Положим
г = -^С(1-С)2.
Для всех четных к = 2т из рекуррентного уравнения путем итераций получаем сумму первых т членов геометрической прогрессии
а2т = 1 + Г + • • • + гт.
В силу того, что знаменатель геометрической прогрессии г удовлетворяет условию
И = - С)2 = - е^с)2 < 1,
подпоследовательность {а2т} сходится при т ^ ж к пределу
12
1 - г 2 + С - 2 С2 + С3 '
Легко проверить, что подпоследовательность |а2т+і| при т ^ ж сходится к тому же самому пределу а. Следовательно, ак ^ а при к ^ ж.
Последовательность функций {Фк} тогда сходится к функции
( 2+0-204^^’ ЄСЛИ 1 ~ °>
Ф(і) = < 2+с-2С2+с3 + С + С3 — 2Се~м*), если 0 < і < с,
I 1, если і > с;
а последовательность функций {Фк} сходится к функции
Ф(і) = - с) =
2 + С-2С2 + С3 (-1 + если 0 < і < с,
2-|-д—2с2-|-с3 С* — (4 + С* + С*3)е ^ -\- 2е ; если с < і < 2с,
1 — в-м*, если і > 2с.
Нетрудно проверить, что Ф и Ф являются функциями распределения некоторых случайных величин, которые будем обозначать соответственно через X и У.
6. Вычисление скорости роста вектора состояний. Сначала заметим, что из сходимости последовательностей функций распределения Фк ^ Ф и Фк ^ Ф следует сходимость Рк ^ р и Цк ^ Ц. Определим функцию
2 Ц.С2
2 + С-2С2 + С3 еесли ^ <1 < С,
ф(г) = Ф’(г) = ( 2+с-2С2+е3 ((4 + С + С3)е-^ - 4е-2^), если с < I < 2с,
если Ь > 2с,
а также функцию
2 цС
Ф(Ь) = Ф'(Ь) =
Вычислим вероятности
р = P{X = с} = 1 — ц = Р{У = 0} = 1 —
2+С-2С2 + С3 2 ЦС 2+С-2С2 + С3 '
0,
—
если Ь < 0, если 0 < Ь < с, если Ь > с.
Ф(Ь)ЙЬ : ф(Ь)йЬ
(2 + С)(1 - С)3 2 + С - 2С2 + С3 ’ (2 + С)(1 - С)2 2 + С - 2С2 + С3 ’
После подстановки С = е мс получим
(2 + е~мс)(1 -е^с)3
2 + - 2е~2^с + е~3^с
Р =
(2 + е^с)(1 -е^с)2 2 + - 2е~2^с + е~3^с
В силу сходимости последовательности функций распределения случайных величин X(к) к функции распределения X имеем сходимость средних EX(к) ^ EX.
Тогда для средней скорости роста вектора состояний системы выполняется
1 1 .к 1 .к
А= Ит -Етах(0, У(/г)) + Ит — У^ЕХ(/г)= Ит — ЕХ(к) = ЕХ.
Ь—1с Ь—к* ■ ^ £•—к* ■ ^
к——о к
Найдем интеграл
С (2 + 6^с — Б^сС + ^сС3) /„ (да<4= М2 + С-2С-+С»)
2С
Теперь можно записать
р О
А = ср + ^ ^(()Д -»+^-+ 6. _ 26, + 6.3) •
Заменяя С на соответствующую экспоненту, окончательно получаем
Л — с +
2е_мс
^(2 + е_^с — 2е_2^с + е_3^с)'
Подставляя с = 0, находим Л = 1/^, что соответствует результату, полученному для этого случая выше.
о
°
°
Рис. 1. Зависимость средней скорости роста Л от величины с и fx.
Графики зависимости величины Л от с и м представлены на рис. 1. В заключение заметим, что при с = 1/м выполняется
А=^,
М
где
2е-1
<?1 = 1 + т;---;---г—5------о « 1-3427.
2 + е-1 - 2е-2 + е-3
Литература
1. Baccelli F., Cohen G., Olsder G. J., Quadrat J.-P. Synchronization and linearity: An algebra
for discrete event systems. Chichester: Wiley, 1992. 514 p.
2. Маслов В. П., Колокольцов В. Н. Идемпотентный анализ и его применение в оптимальном управлении. М.: Физматлит, 1994. 144 с.
3. Heidergott B., Olsder G. J., van der Woude J. Max-Plus at work: Modeling and analysis of synchronized systems. Princeton: Princeton University Press, 2006. 226 p.
4. Кривулин Н. К. Методы идемпотентной алгебры в задачах моделирования и анализа
сложных систем. СПб.: Изд-во С.-Петерб. ун-та, 2009. 256 с.
5. Olsder G. J., Resing J. A. C., De Vries R. E., Keane M. S., Hooghiemstra G. Discrete event systems with stochastic processing times // IEEE Trans. Automat. Contr. 1990. Vol. 35, №3. P. 299-302.
6. Jean-Marie A. Analytical computation of Lyapunov exponents in stochastic event graphs // Performance Evaluation of Parallel and Distributed Systems. Solution Methods: Proc. 3rd QMIPS Workshop. Amsterdam: CWI, 1994. P. 309-341. (CWI Tracts, Vol. 106.)
7. Кривулин Н. К. Вычисление показателя Ляпунова для одной модели обобщенной линейной стохастической динамической системы // Математические модели. Теория и приложения. Вып. 10. Сб. науч. статей / Под ред. М. К. Чиркова. СПб.: ВВМ, 2009. С. 105-110.
8. Кривулин Н. К. Вычисление показателя Ляпунова для одной модели стохастической системы с синхронизацией // Стохастическая оптимизация в информатике. Вып. 6. / Под ред. О. Н. Граничина. СПб.: Изд-во С.-Петерб. ун-та, 2010.
Статья поступила в редакцию 7 сентября 2010 г.