Статистическое моделирование гауссовской очереди Текст научной статьи по специальности «Математика»

Труды Карельского научного центра РАН № 5. 2011. С. 55-62

УДК 519.872.6

СТАТИСТИЧЕСКОЕ МОДЕЛИРОВАНИЕ ГАУССОВСКОЙ ОЧЕРЕДИ

О. В. Лукашенко1, Е. В. Морозов1, М. Пагано2

1 Институт прикладных математических исследований Карельского научного центра РАН

2 Факультет информационного инжиниринга,

Университет Пизы, Италия

Статья посвящена оцениванию некоторых характеристик гауссовских очередей (систем с гауссовским входным потоком). Основное внимание уделено системе с фрактальным броуновским входным процессом (фрактальное броуновское движение, ФБД). Исследуется аппроксимация вероятности большого уклонения процесса загрузки. Методом имитационного моделирования анализируется гауссовская очередь в режиме большого числа входящих потоков. Работа выполнена при поддержке РФФИ, проект 10-07-00017.

Ключевые слова: гауссовская очередь, структура зависимости, фрактальное броуновское движение.

O. V. Lukashenko, E. V. Morozov, M. Pagano. SIMULATION OF THE GAUSSIAN QUEUE

The paper is devoted to the dependence structure of Gaussian queues (queues fed by Gaussian input). The focus is on the system with fractional Brownian input (fractional Brownian motion, fBm). An approximation of the probability of large deviation is investigated. A Gaussian queue in the multiple-source mode is analyzed by simulation. The study was supported by RFBR, project 10-07-00017.

Key words: Gaussian queue, dependence structure, fractional Brownian motion.

Введение

В настоящее время гауссовские модели систем обслуживания достаточно широко используются, поскольку позволяют адекватно описать входной поток для широкого класса сетевых узлов ^опте, 1994; КПр1, Коггов, 2002; Mandjes, 2007]. В частности, они позволяют включить в моделирование самоподобие и долговременную зависимость, которые играют важную роль в описании динамики современных телекоммуникационных систем. Напомним, что самоподобие означает, что по-

ведение процесса качественно не меняется при изменении масштаба времени, а долговременная зависимость означает столь медленное убывание автокорреляционной функции, что автокорреляционный ряд оказывается несум-мируемым. Самым известным и наиболее хорошо изученным гауссовским процессом является фрактальное броуновское движение (ФБД), которое обладает свойствами самоподобия и долговременной зависимости. В частности, ФБД возникает при суперпозиции большого числа независимых так называемых оп/ой-источников с тяжелыми хвостами на

■©

больших масштабах времени. (Более подробно см. [Taqqu et al., 1995].)

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

P (Q > b) ~ fi(b), b ^ те P (Qn > nb) ~ f2(n), n ^ те,

где a ~ b означает, что a/b ^ 1, называются точными асимптотиками, здесь fi, f2 - некоторые явно заданные функции, Q - процесс загрузки. Иногда удается получить лишь логарифмические асимптотики:

ln P (Q > b) ~ fe(b), b ^ те ln P (Qn > nb) ~ f4(n), n ^ те,

дающие менее полную информацию о вероятности переполнения. Отметим наиболее важные работы [Duffield, O’Connell, 1995; Narayan, 1998; Husler, Piterbarg, 1999; Debicki, Mandjes, 2002], посвященные этой проблематике. В работе [Es-Saghouani, Mandjes, 2009] предпринята попытка исследовать структуру зависимости гауссовской очереди. Поскольку выражение для корреляционной функции гауссовской очереди в явном виде получить пока не удалось, в [Es-Saghouani, Mandjes, 2009] введена и исследована иная (но близкая) мера зависимости, также основанная на двумерном распределении процесса загрузки. Для этой функции найдена асимптотика для случая большого числа источников. В данной статье мы постарались найти такое количество источников (приемлемый уровень агрегирования), когда указанные асимптотические результаты имеют достаточную точность. В статье также приведены результаты имитационного моделиро-

вания по оценке вероятности превышения заданного уровня в гауссовской очереди с ФБД в качестве входного потока.

ГАУССОВСКИЕ ОЧЕРЕДИ

Пусть A(t) есть объем трафика, поступающего в систему на интервале времени [0,t], и пусть A(s,t) = A(t) — A(s), t > s > 0. Заметим, что A(—t) =d — A(t), где =d означает равенство по распределению. Будем предполагать, что входной поток (входной трафик) имеет следующий вид [Norros, 1994]:

A(t) = mt + л/атБн (t), (1)

где m > 0 - средняя интенсивность трафика, Bh(t) - ФБД, описывающее флуктуации трафика вокруг линейно возрастающего среднего, 1/2 < H < 1 - так называемый параметр Херста, а > 0 - некоторая константа. Предполагается, что система имеет одно обслуживающее устройство с постоянной скоростью обслуживания C > т. Процесс загрузки (в стационарном режиме) для такой системы определяется следующим образом (см., например, [Mandjes, 2007]):

Qt =d sup (A(s, t) — C(t — s)),

s^t

Ввиду стационарности Qt =d Qs для любых t, s. Для удобства иногда зависимость Qt от t будем опускать. Далее предполагаем, что имеется n независимых, одинаково распределенных гауссовских источников описанного выше вида, а скорость обслуживания пропорциональна количеству источников, т. е. C = nc. (Это обычное предположение при асимптотическом анализе системы с растущим числом источников.) Тогда с учетом (1) стационарный процесс загрузки в момент времени t удовлетворяет следующим стохастическим равенствам:

Aj(s, t) — nc(t — s)^j

(n

i/amBH(s,t) — nr(t — s)

1=1

=d sup ( \/amnBH (s,t) — nr(t — s))

s^t

=d sup (i/amnBH (s) — nrs) ,

s^0

здесь BH(u, v) = BH(v)—BH(u) - нагрузка, поступающая от источника i на интервале [u, v], u ^ v, B(u,v) - нагрузка от типичного источника, а параметр r = c — m. Отметим, что

Qtn =d sup

s<t

предполагается выполненным условие г > 0, обеспечивающее существование стационарного режима для процесса загрузки.

Для произвольного заданного Ь стационарная вероятность переполнения (превышения уровня В = пЬ) определяется как

nn(b) = P (Qn > nb).

(2)

І

\/2п

e-x2/2dx. (4)

В общем случае (1/2 < H < 1) аналогичный точный результат пока не найден. Однако имеет место ряд асимптотических результатов. В частности, было доказано [Botvich, Duffield, 1995], что справедливы следующие предельные соотношения:

1 V 2(t)

— lim — ln nn(b) = inf ---------, (5)

n^ro n teK+ 2am

где

V(t) =

(b + rt)2 tH

Значение £*, при котором достигается минимум в (5) - это так называемое наиболее вероятное время переполнения. Для ФБД можно получить, что

Точные выражения для (2) получены лишь в некоторых простейших случаях. Так, при Н = 1/2, т. е. для стандартного броуновского движения, в работе [Такач, 1971] доказано, что

Пп(Ь) = ехр Г -22ПГ^^ .

\ ат у

Также точное выражение можно получить и при Н = 1. В этом случае процесс Ві(і) = і ■ N(0,1) является случайной прямой (лучом в случайном направлении из начала координат), где N(0,1) - стандартная нормальная случайная величина. В этом случае имеем

Qn =d sup ( \/amnB1(s) — nrs)

s^0

=d sup ((\/amnN(0,1) — nr) s) .

s^0

Из этого следует, что

n [те, если N(0,1) > ,

Qn = < ’ v ’ ' (3)

0, иначе

Выражение (3) имеет понятную физическую интерпретацию. В систему в каждый момент времени поступает одинаковая нагрузка, поэтому возможны два случая: либо система справляется с этой работой, а значит, очереди не образуется, либо работа неограниченно накапливается. Тогда для вероятности переполнения nn(b) при любых b справедливо следующее выражение:

nn(b) = P (n(0,1) > //£)

am

t=

Hb

V(t*) =

І - H r

1H

(б)

b

1 - H

r \H

~H)

Выражение (5) означает, что при достаточно большом n

nn(b) ~ exp (-2amv2(t,^. (7)

Структура зависимости процесса загрузки

Для изучения структуры зависимости проанализируем двумерное распределение процесса загрузки. Естественный подход к изучению такой структуры - это анализ корреляционной функции:

<№) =

EQnQn - (EQn)(EQn) ^Var(Qn)Var(Qn) '

Очевидно, что 5га(Т) ^ 0 при Т ^ то. Однако остается неясной скорость убывания, так как для ФБД точная формула даже для средней стационарной нагрузки не найдена.

Поэтому вводятся в рассмотрение иные меры зависимости. Мы будем анализировать следующую характеристику меры зависимости, введенную в работе [Es-Sаghouаni, Мапфев, 2009]:

Kn(T) =

P(Qn > np, Qn > nq) P(Qn > np)P(Qn > nq):

(8

где параметры р, д > 0. В [Es-Saghouani, Mandjes, 2009] также показано, что для (Т) справедлива асимптотика:

lim — in Kn(T) = kGT n——<^ n

2H-2

+ o(T2H-2), (9)

где параметры

kG =

p

(2H — і)г2 s2-2H+2-

Hm2 *

2H

H

s* —

1- H

t* —

H

1 - H’

r

r

(см. (6)). Это означает, что при достаточно больших n

Kn(T) и exp (nkoT2H-2 + n ■ o (T2H-2)) .

(10)

На первый взгляд, эта аппроксимация может быть использована для оценки зависимости процесса загрузки. Однако формулу (10) нужно применять с осторожностью, ввиду присутствия множителя n при неизвестной функции o (T2Н-2).

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

Л(Т ) =

P (WGn > np, Wy > nq)

P (W0n > np) P (Wy > nq),

где

W0n = £ Aг(0,1), W? = £ Aг(T,T + 1),

г=1

г=1

т. е. суммарная работа на интервалах [0,1] и [Т, Т + 1] соответственно, причем в силу стационарности приращений ФБД, = ^р. Обозначим

Л(Т) = 11ш — 1п Л(Т).

П

По теореме Крамера (см. [Mandjes, 2007; Ев-Saghouani, Mandjes, 2009])

Это позволяет говорить о том, что долговременная зависимость входного потока ФБД влечет за собой долговременную зависимость процесса загрузки.

Оценивание стационарных характеристик

В данном разделе рассматриваются оценки величин nn(b),Kn(T). Для этого генерируется N реализаций ФБД, каждая из которых содержит M наблюдений (т. е. M -длина траектории). Для построения этих реализаций использован RMD-метод (Random Midpoint Displacement) [Lau et al., 1995], временная сложность которого составляет O(M). Для оценивания стационарных характеристик системы, длина M генерируемой траектории должна быть достаточно большой, чтобы устранить (или по крайней мере уменьшить) влияние начального, так называемого переходного периода. (Отметим, что подробно вопрос моделирования ФБД обсуждается в [Coeurjolly, 2000].)

Обозначим k-ю последовательность (траекторию) значений ФБД X

1, ... , N. Обозначим также

XjM, где к =

л/anmXk(s,u) — nr(u — s) , (12

где £ = {0,1,...,и}. Стандартная выборочная оценка для пр(Ь) имеет вид:

Л(Т) ~ const ■ T2Н 2, T ^ те,

что согласуется с (9).

Также можно говорить о наличии долговременной зависимости (долгой памяти) для процесса загрузки при достаточно больших n и T. Действительно, обозначим

Rn(T) = Corr (/(Qn > np), I(Qn > nq)),

где / означает индикатор. Принимая во внимание соотношения (7), (10), можно получить (см. [Mandjes, 2007; Es-Saghouani, Mandjes, 2009]), что при достаточно больших значениях параметров n и T

Rn(T) и (enkoT2H-2 — 1)^(n) и ^(n)T2H-2,

(11)

где ^(п) - некоторая (независящая от Т) функция. Выражение (11) показывает, что функция Яп(Т) несуммируема при Н > 1/2.

N

£i(Q U > nb

г=1

где момент u выбран достаточно большим.

Sample size

Рис. 1. Стабилизация оценки Пю(1) с ростом числа наблюдений

n

Рис. 2. Сравнение двух методов оценки стационарной вероятности переполнения пп (1)

n

Рис. 3. Зависимость вероятности переполнения пп(1) от количества источников

Рис. 1 показывает зависимость оценки вероятности переполнения піо (1) от длины траектории. При моделировании использованы следующие значения параметров: п = 10; Н = 0, 8; с = 1; т = 0, 8; Ь = 1; N = 104. В соответствии с формулой (6) получаем і* = 20. На рис. 1 видна определенная стабилизация значений оценки после наиболее вероятного времени переполнения і*. Отметим вновь, что для оценивания стационарных характеристик начальная часть траектории должна быть пропущена, т. е. общая длина траектории должна быть достаточно большой.

Стационарную вероятность переполнения (2) можно оценивать и как долю времени нахождения процесса загрузки выше уровня В = пЬ, т. е.

M 21 (Q(kh) >nb), (13)

k=1

где процесс загрузки в дискретном времени определяется известной рекурсией Линдли:

Q(t + h) = (Q(t) — nhr + /namBH (t, t + h))+ ,

где h > 0 - так называемый шаг дискретизации (выбор h обсуждается ниже), Q(0) = 0.

Значение (13) вычисляется для каждой траектории. В качестве оценки вероятности переполнения берется среднее полученных значений. Заметим, что совпадение в пределе этих оценок теоретически доказано лишь на основе так называемого свойства PASTA, которое применимо лишь в случае пуассоновско-го входного потока. Однако рис. 2 показывает, что результаты, полученные по схеме (12) (Simulation1) и схеме (13) (Simulation2), не отличаются значительно и для рассматриваемой модели.

Следует отметить, что схема (13) предпочтительнее с точки зрения вычислительных затрат. Так, например, используя N = 1000 траекторий ФБД, этот метод позволяет оценить вероятности порядка 10-5 (см. рис. 3), в то время как метод (12) позволяет достичь такой же точности при генерации как минимум N = 105 траекторий, что естественным образом сказывается на времени работы алгоритма.

Теперь сравним результаты численного моделирования вероятности nn(b) с аппроксимацией (7). Чтобы определить, при каком n эта аппроксимация имеет заданную точность, рассмотрим так называемые невязки

An = |nn(b) — Пп(Ь)| , (14)

при некоторых фиксированных значениях Ь, г. Требуется найти такое минимальное значение п, при котором Др < е, для некоторого заранее заданного значения е > 0. (Таким образом, е - это точность аппроксимации.)

На рис. 4 представлен график зависимости Дп от числа источников п (логарифмический масштаб) при следующих параметрах: Ь = 1; с = 1; т = 0, 8; Н = 1; N = 1000; М = 10000. Этот график наглядно показывает, каким должен быть уровень агрегирования п, чтобы обеспечить необходимую точность аппроксимации (7).

ный интерполяционный сплайн, построенный по точкам Xo, Xh,..., XMh. В этой связи уместно привести аналитический результат о величине ошибки такой дискретизации в случае стандартного броуновского движения B^2(t).

Пусть h = hN = 1/L, где L > 0, Bh(t) - реализация стандартного броуновского движения с шагом h. Тогда средняя ошибка линейной интерполяции удовлетворяет соотношению (см. [Asmussen, Glynn, 2007]):

Рис. 4- Зависимость Дп от числа источников Е Bh(t) — Bi/2(t) dt = л/n/32/L1/2.

Jo

Зависимость Дп от шага дискретизации

№ Дп, h = 1 СО 1 1 h= ,n <

1 0,151152 0,150463

2 0,148348 0,149770

3 0,148442 0,146604

4 0,150842 0,146778

5 0,151446 0,151623

6 0,151327 0,147263

7 0,148494 0,144033

8 0,152131 0,147650

9 0,147023 0,148453

10 0,151670 0,149856

11 0,148143 0,150300

12 0,149802 0,145306

13 0,149432 0,148748

14 0,155012 0,145321

15 0,147368 0,148984

16 0,152423 0,151441

17 0,147264 0,149285

18 0,148313 0,150049

19 0,147407 0,150392

20 0,150885 0,152498

Среднее: 0,149846 0,148741

Результаты эксперимента представлены в таблице, где видна тенденция к уменьшению Дп с уменьшением Н. Таким образом, чтобы точнее оценивать характеристики случайного процесса с непрерывным временем, необходимо уменьшать шаг дискретизации Н. С другой стороны, уменьшение Н ведет к увеличению объема вычислений, что ставит проблему баланса между величиной шага Н и объемом вычислений.

На рис. 5 представлен график зависимости (логарифмический масштаб) оценки вероятности переполнения от величины параметра Херста Н при следующих значениях остальных параметров: Ь = 1; с = 1; т = 0, 8; п = 10; М = 104; N = 104.

На рис. 6 представлен график зависимости вероятности переполнения от размера буфера Ь при Н = 1 со следующим параметрами: п = 10; с = 1; т = 0, 8; М = 104; N = 104. В силу (4) П1о(Ь) 0, 23975 для любого Ь. Как

видно, результаты численного моделирования согласуются с теоретическими.

Для изучения влияния шага дискретизации Н на результаты оценивания, для двух значений Н = 1/16, 1 проведена серия из 20 испытаний. Каждое испытание состоит в подсчете значения выражения (14) по схеме (13) при фиксированных п = 10 и Ь = 1. Поскольку аппроксимация (7) получена для процесса загрузки в непрерывном времени, то можно ожидать, что с уменьшением величины шага Н значения Дп также должны уменьшаться. Отметим, что при численном моделировании исходный процесс с непрерывным временем заменяется на процесс с дискретным временем, что неизбежно приводит к ошибке дискретизации. Фактически траектория процесса представляет собой линей-

Рис. 5. Зависимость Пю(1) от параметра Херста H

n

H

0

10 15 20

Buffer size b

Рис. б. Оценка вероятности nio(b) при H =1

Чтобы оценить величину КП(Т), нужно оценить вероятности, входящие в выражение (8). Напомним, что величины <50, 5т определяются в (12). Тогда оценка двумерного распределения Р (5П > пр, > пд) определяется как выборочное среднее:

1 ^

Пп,т(р, 5) = N ^ I к > пр, > п<?) . к=1

С учетом этого оценка для Кп(Т) запишется в виде отношения:

Kn(T) =

nn,T (P,g) 7Tn(p)7rn(g) '

Рис. 7. Оценка K10(T)

На рис. 7 представлен результат численного моделирования величины Кю (Т) и его сравнение с аппроксимацией (10). В моделировании использованы следующие значения параметров: Н = 0, 8; с = 1; т = 0, 8; п = 10; р = 1; д = 2; N = 104.

Поскольку получение большого количества траекторий ФБД большой длины за приемле-

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

Заключение

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

Литература

Такач Л. Комбинаторные методы в теории случайных процессов. М.: Мир, 1971. 264 с.

Asmussen S., Glynn P. Stochastic Simulation: algorithms and analysis. Springer. 2007. P. 488.

Botvich D., Duffield N. Large deviations, the shape of the loss curve, and economies of scale in large multiplexers // Queueing Systems. 1995. Vol. 20. P. 293-320.

Coeurjolly J. F. Simulation and identification of the fractional Brownian motion: a bibliographical and comparative study // Journal of Stat. Software. 2000. Vol. 5.

Debicki K., Mandjes M. Exact overflow asymtotics for queues with many Gaussian inputs. Report PNA-R0209 March 31, 2002.

Duffield N., O’Connell N. Large deviations and overflow probabilities for the general single server queue, with applications // Proceedings of the Cambridge Philosophical Society. 1995. Vol. 118. P. 363-374.

Es-Saghouani A., Mandjes M. On the dependence structure of Gaussian queues // Stochastic Models. 2009. Vol. 25. P. 221-247.

Husler J., Piterbarg V. I. Extremes of a certain class of Gaussian processes // Stochastic Processes and their Applications. 1999. Vol. 83. P. 257-271.

Kilpi J., Norros I. Testing the Gaussian approximation of aggregate traffic, Proceedings of the second ACM SIGCOMM Workshop, Marseille, France. 2002. P. 49-61.

Lau W-C., Erramilli A., Wang J. L., Willinger W. Self-similar traffic generation: the random midpoint displacement algorithm and its properties. Proceedings of ICC ’95. 1995.

Mandjes M. Large Deviations of Gaussian Queues. Chichester: Wiley, 2007. P. 339.

Narayan O. Exact asymptotic queue length distribution for fractional Brownian traffic // Advances in Performance Analysis. 1998. Vol. 1. P. 39-63.

O

O

5

25

30

T

б1

Norros I. A storage model with self-similar input a fundamental result in self-similar traffic modeling

// Queueing Systems. 1994. Vol. 16. P. 387-396. // Computer communication review. 1997. Vol. 27.

Taqqu M. S., Willinger W, Sherman R. Proof of P. 5-23.

СВЕДЕНИЯ ОБ АВТОРАХ:

Лукашенко Олег Викторович

аспирант

Институт прикладных математических исследований КарНЦ РАН

ул. Пушкинская, 11, Петрозаводск,

Республика Карелия, Россия, 185910 эл. почта: lukashenko-oleg@mail.ru тел.: (8142) 763370

Морозов Евсей Викторович

ведущий научный сотрудник д. ф.-м. н.

Институт прикладных математических исследований КарНЦ РАН

ул. Пушкинская, 11, Петрозаводск,

Республика Карелия, Россия, 185910 эл. почта: emorozov@karelia.ru тел.: (8142) 763370

Пагано Микеле

преподаватель

факультет информационной инженерии Университета

г. Пизы, Италия

эл. почта: m.pagano@iet.unipi.it

тел.: +39 050 2217575

Lukashenko, Oleg

Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Science 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: lukashenko-oleg@mail.ru tel.: (8142) 763370

Morozov, Evsey

Institute of Applied Mathematical Research, Karelian Research Centre, Russian Academy of Science 11 Pushkinskaya St., 185910 Petrozavodsk, Karelia, Russia

e-mail: emorozov@karelia.ru tel.: (8142) 763370

Pagano, Michele

Associated Professor

Department of Information Engineering,

University of Pisa, Italy e-mail: m.pagano@iet.unipi.it tel.: +39 050 2217575