Особенности моделирования систем массового обслуживания с тяжелыми хвостами распределений на GPSS World. Метод Arand Текст научной статьи по специальности «Математика»

ФИЗИКО-МАТЕМАТИЧЕСКИЕ НАУКИ

уДК 5192:004.421.5 004.7 в. Н. ЗАДОРОЖНЫЙ

Омский государственный технический университет

ОСОБЕННОСТИ МОДЕЛИРОВАНИЯ СИСТЕМ МАССОВОГО ОБСЛУЖИВАНИЯ С ТЯЖЕЛЫМИ ХВОСТАМИ РАСПРЕДЕЛЕНИЙ НА GPSS WORLD. МЕТОД ARAND

Рассматриваются проблемы построения адекватных математических моделей обслуживания фрактального трафика, актуальные для современной теории проектирования компьютерных сетей. Для решения проблемы корректной реализации распределений с тяжелыми хвостами предлагается метод ARAND (Accurate RAND). Выявляются дефекты механизма продвижения времени в GPSS World. Разрабатываются рекомендации по моделированию на GPSS фрактальных систем массового обслуживания.

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

1. Введение. Имитационное моделирование — дефекты механизма продвижения времени

(ИМ) на GPSS World [1] систем массового обслужи- в GPSS, приводящие к некорректности оценок, по-

вания (СМО), характеризуемых распределениями лучаемых при длинных прогонах;

с тяжелыми хвостами (РТХ), затрудняют следую- — несовместимость разных ГСЧ при длинных

щие неочевидные, но серьезные проблемы: прогонах;

— короткий период имеющихся в GPSS генера- — искажение РТХ при их реализации в ИМ [2]. торов стандартных случайных чисел (ГСЧ), недоста- Первые три проблемы всплывают и при моточный для реализации длинных прогонов; делировании на GPSS «классических» СМО,

характеризуемых распределениями с легкими хвостами (РЛХ), если нужны длинные прогоны [3]. Проблема искажения РТХ характерна и для других систем ИМ, хотя проявляется в них менее заметно благодаря более высокой, чем в GPSS, разрядности используемых ГСЧ.

В статье приводятся результаты исследования вышеуказанных проблем и разрабатывается методика ИМ СМО с РТХ на GPSS World, обеспечивающая высокую точность получаемых оценок. Предлагается ГСЧ ARAND (Accurate RAND), обеспечивающий корректную реализацию РТХ в любых системах ИМ.

2. О длине периода генераторов случайных чисел GPSS. Встроенные в GPSS World ГСЧ, согласно [1], основаны «на мультипликативно-конгруэнтном алгоритме Лемера с максимальным периодом. Алгоритм генерирует псевдослучайные числа в открытом интервале от 0 до 2 147 483 647 и до самоповторения генерирует 231 — 2 = 2 147 483 646 уникальных псевдослучайных чисел. Дополнительно в GPSS World используется шаг перемешивания». Однако несложные, хотя и трудоемкие эксперименты, выполненные на GPSS World, показали следующее.

1. Последовательность псевдослучайных чисел на выходе процедуры Uniform(1,0,1) — первого ГСЧ в GPSS World — начинает самоповторяться с числа номер 2 048 000 000 = 211х106, т.е. длина его периода равна 2 047 999 999 = 8 863x743x311. Эти числа, как и указано в [1], принимают значения от 0,000000 до 0,999999 (все разряды дробной части после шестого равны нулю).

2. Периоды длиной 2 047 999 999 чисел имеют также ГСЧ Uniform(2,0,1), Uniform(3,0,1) и Uniform(4,0,1). Эти ГСЧ выдают последовательность чисел из одно-

го и того же направленного кольца (при этом номер j ГСЧ Uniform(j,0,1) определяет, какой элемент кольца становится началом последовательности [1]).

Таким образом, сведения о ГСЧ GPSS World, приведенные в [1], не вполне соответствуют действительности. Но порядок длины периода ГСЧ указан верно.

Знание длины периода используемых ГСЧ необходимо для корректной организации длинных прогонов имитационных моделей.

3. Дефекты времени в GPSS World. В статье [3] приведены результаты моделирования на GPSS системы M/M/1 с интенсивностью входящего потока X = 1 и коэффициентом загрузки р = 0,9. Средняя длина очереди L = 8,1 и среднее время ожидания w = 8,1 для этой системы известны точно. Результаты ИМ этой СМО при различной длине прогона, выраженной (средним) числом N заявок, поступающих за моделируемое время TM, представлены в табл. 1.

Во второй строке таблицы приведены абсолютные погрешности полученных имитационных оценок w времени w. В третьей (последней) строке приведена допустимая погрешность, соответствующая правилу трех сигм: она найдена методом, предложенным в [3] для выборок с коррелирующими элементами w.. В GPSS-модели интервалы поступления и обслуживания заявок заданы блоками GENERATE (Exponential(1,0,1)) и ADVANCE (Exponential(1,0,0.9)). Реализующая экспоненциальную случайную величину (с.в.) функция Exponential (j,0,E) использует очередное число, выдаваемое ГСЧ Uniform(j,0,1) (обозначим его через z), преобразует его по формуле x = — Eln(z) (попадающиеся нулевые значения z пропускаются), и возвращает x.

Результаты моделирования системы M/M/1 на GPSS World

Таблица 1

Показатель Число испытаний N = VT., M

1 млн 10 млн 100 млн 1 млрд 2 млрд

Оценка среднего времени ожидания 8,332 8,061 8,087 8,043 7,799

Фактическая погрешность оценки 0,232 0,039 0,013 0,057 0,301

Допустимая погрешность («три сигмы») 0,404 0,064 0,040 0,006 0,004

GENERATE (Exponential(1,0,1))

ASSIGN 1,AC1

QUEUE 1

SEIZE 1

DEPART 1

SAVEVALUE 1+,((AC1-P1)#(AC1<P1))

SAVEVALUE 2+,(AC1<P1)

ADVANCE (Exponential(1,0,0.9))

RELEASE 1

TERMINATE

GENERATE 100000000

SAVEVALUE W NEG,(X1/X2)

TERMINATE 1

START 1

RESET

START 1

RESET

START 1

Таблица 2

Характеристики отрицательных значений времени ожидания

Рис. 1. Программа, выявляющая наличие обратного хода времени в GPSS World

Абсолютное время AC1 Число случаев wt < 0 Сумма всех wt < 0 Среднее значение отрицательных w

1x10" 344 -9,87597 -0,028709

2x10" 4961 -294,378 -0,059338

3x108 23470 -2085,86 -0,088873

4x108 69841 -8161,87 -0,116864

5x108 161040 -23295,5 -0,144657

6x108 317228 -54594,3 -0,172098

7x108 557714 -110980 -0,198991

8x108 901127 -202575 -0,224802

9x108 1354436 -337396 -0,249104

109 1920013 -522135 -0,271943

GENERATE (Exponential 1 0 1))

QUEUE 1

SEIZE 1

DEPART 1

ADVANCE (Exponential 1 0 0.9))

RELEASE 1

TERMINATE

GENERATE 100000000

TERMINATE 1

START 1

CLEAR

START 1

CLEAR

START 1

Рис. 2. Метод увеличения объема выборки при ИМ СМО с РЛХ

Описанный M/M/1-тест показал, что, независимо от номера j применяемого ГСЧ, при длине прогона N > 100 + 200 млн результаты ИМ некорректны (табл. 1), в соответствии с чем в [3] дана рекомендация такие прогоны не использовать.

Однако настоятельная необходимость длинных прогонов при ИМ СМО с РТХ требует выяснения причин смещения оценок, возрастающего с ростом длины прогона (табл. 1). Выполненные для этого эксперименты (описанные ниже) привели к открытию малоприятного факта: механизм продвижения времени в GPSS World имеет скрытые дефекты.

На рис. 1 представлена программа, выполнение которой доказывает наличие дефектов времени в GPSS World. В ходе выполнения этой программы, моделирующей СМО M/M/1, ячейка 1 накапливает сумму всех отрицательных значений w . времени ожидания транзактов. В ячейке 2 подсчитывается число отрицательных w.. Среднее значение отрицательных w. вычисляется в ячейке W_NEG. Значения ячеек выводятся после каждого из 10 отрезков длинного прогона, запускаемых командами START. В каждом отрезке прогона модельное время AC1 продвигается вперед на 100 млн единиц времени.

Результаты выполнения этой программы представлены в табл. 2.

Незначительное удлинение программы позволяет выяснить, что первый транзакт с отрицательным временем ожидания — это транзакт с номером 21 693 107. Он входит в очередь в момент модельного времени 21 693 041,901 257, затем выходит из нее в момент 21 693 041,894 325, который меньше времени входа на Д1 = 0,006 932 (т.е. время «идет вспять»). Заметим, что в среднем отрицательные w. растут (по модулю) пропорционально времени AC1 (табл. 2), а частота появления отрицательных w и модуль их суммы растут примерно пропорционально четвертой-пятой степени времени AC1. Наличие отрицательных w доказывает некорректность работы планировщика транзактов со списками событий. В целом же искажаются и положительные w.. Эксперименты показывают, что вместе с длиной прогона растут искажения и других оценок вероятностно-временных характеристик (ВВХ): средней длины очереди, вероятности потери заявки (при ограниченной длине очереди) и т.д.

С внешними проявлениями описанного дефекта разработчики GPSS World, по-видимому, сталкивались, но сам дефект не распознали. Мы можем судить об этом по тому, например, как в [1] (с. 306)

смутно замечается, что использование команды RESET «может привести к некоторой систематической ошибке в области наименьших значений». Однако табл. 2 ясно показывает, что первые 344 отрицательных w появляются в нашем тесте задолго до выполнения первой команды RESET, так что она здесь ни при чем. Проверки показывают, что обратный ход времени не обусловлен также и погрешностями округления.

Экспериментально мы установили, что в AnyLogic [4] такого дефекта нет.

4. Реализация длинных прогонов при моделировании СМО с РЛХ. Объем выборки при моделировании СМО с РЛХ на GPSS World можно все же увеличить до пропускания через модель 1 млрд тран-зактов (т.е. до 2 млрд обращений к ГСЧ Uniform), не опасаясь ущерба от дефектов модельного времени. Для этого длинный прогон достаточно разбить командами CLEAR и START на 10 независимых прогонов умеренной длины, составляющей 100 млн транзактов, как это показано в варианте M/M/1-теста на рис. 2.

При этом дефекты модельного времени не успевают причинить заметный вред, поскольку абсолютное модельное время «своевременно» обнуляется командами CLEAR. Кроме того, благодаря быстрому затуханию переходных процессов (ПП) в СМО с РЛХ, все прогоны умеренной длины дают корректные оценки стационарных ВВХ. А так как эти оценки независимы, то, усредняя их, нетрудно оценивать точность получаемого при усреднении результата.

Например, погрешность получаемой таким образом итоговой оценки для w в СМО M/M/1 при р = 0,9 (см. рис. 2) находится в допустимых пределах. Отдельные 10 прогонов длины 100 млн единиц времени дают для w оценки \v1,...,w10 = 8,120, 8,070, 8,070, 8,111, 8,078, 8,134, 8,128, 8,113, 8,087, 8,118. Каждая оценка отклоняется от точного значения 8,1 менее чем на максимально допускаемую правилом трех сигм величину 0,04 (для прогонов такой длины). Накапливаемые по прогонам средние Wk = (1/kwt , — это десять оценок Я,..., wi0 = 8,120, 8,095, 8,087, 8,093, 8,090, 8,097, 8,102, 8,103, 8,101, 8,103. Их отклонения от 8,1 также не превышают порогов, определяемых для них как 0,04 /4k . Метод можно применять при любом номере j ГСЧ (на рис. 2 j = 1).

В прогонах умеренной длины дефекты модельного времени на точности оценок стационарных ВВХ еще не сказываются. И уже не сказывается возобновление ПП после каждого CLEAR, поскольку в СМО с РЛХ ПП весьма непродолжительны. Однако для СМО с РТХ этот метод может потребовать доработки, ведь в них даже при небольших р ПП могут длиться очень долго [5].

5. Несовместимость ГСЧ GPSS World в длинных прогонах. Можно попытаться увеличить длину прогона вдвое (до 2 млрд транзактов), затрачивая на каждый транзакт не два обращения к первому Uniform (по одному в блоках GENERATE и ADVANCE, см. рис. 2), а только одно, например, в блоке GENERATE; а в блоке ADVANCE обращаться ко второму Uniform.

На самом же деле, используя в программе на рис. 2 блоки GENERATE (Exponential(1,0,1)) и ADVANCE (Exponential(2,0,0.9)), мы, по причине некой интерференции последовательностей первого и второго Uniform, уже в первом прогоне (выполненном по первой команде START), получим для w неприемлемую оценку 7,898. Ее погрешность

|7,898 — 8,1| = 0,202 более чем впятеро превышает максимально допустимую 0,04 (см. табл. 1 и рис. 2).

Множество дополнительных экспериментов показало, что в длинных прогонах любые ГСЧ GPSS World с разными номерами в общем случае несовместимы. При этом величина искажения результатов ИМ СМО из-за интерференции последовательностей разных ГСЧ зависит от средней длины очереди, от вероятностных распределений, задающих СМО, от ряда других факторов и в общем случае непредсказуема.

6. Реализация длинных прогонов при моделировании СМО с РТХ. При моделировании СМО с РТХ наиболее подходящей организацией имитационного эксперимента на GPSS становится реализация r независимых прогонов длины N < 100 млн заявок. Учитывая, что на каждый транзакт приходится два обращения к ГСЧ, а общее число обращений не должно превышать длины периода ГСЧ 2047999999 » 2 048 млн, имеем ограничение:

r-N < 1,024х109.

(1)

Если ПП вычисляемой оценки (например, и) заканчивается достаточно быстро, то все ее реализации (например, иИ1,...,иг ), полученные в разных прогонах, затем усредняются, причем независимость этих реализаций позволяет обычным образом строить доверительный интервал для итоговой оценки. Если ПП вычисляемой оценки к концу прогонов не заканчиваются, то можно усреднить все г полученных ПП и найти по усредненному ПП оценку стационарного значения показателя методом степенного прогноза [5].

Так, при ИМ СМО Ра/М/1 с параметрами Па-рето-распределения К = 1, а = 1,1 (т.е. со средним временем 11 между приходами заявок) и при среднем времени обслуживания 4,4 (т.е. при коэффициенте загрузки р = 4,4/11 = 0,4) сформулированные выше рекомендации были реализованы следующим образом. Выполнено г = 901 прогонов длины Т = 110 000 единиц времени. Всего приходит в среднем N = г(Т/11) = 901x110 000/11 = 9,01 млн заявок, и ограничение (1) выполняется. В ходе каждого прогона значения ПП оценки и (т.е. системного числового атрибута ОТ1) записывались в 2000 точках отсчета I. = 55, 110, ..., 110000, равномерно расставленных в модельном времени t. Результат ОТ1(Ц усреднения всех 901 ПП графически представлен на рис. 3.

Как видно из рис. 3, ПП здесь заканчивается достаточно быстро, и для оценки стационарного ш можно просто взять последнее значение усредненного ОТ1, составившее здесь 118,4. Этот результат

отличается от точного значения w = 118,5, найденного численным методом, описанным в [5], лишь на 0,07 %.

Такая хорошая точность оценок, получаемых при рекомендуемой здесь организации ИМ, обеспечивается, в частности, тем, что при использовании независимых прогонов отсутствует положительная корреляция между усредняемыми величинами (в силу их независимости). Заметим, что при ИМ этой же СМО (на AnyLogic), в прогоне с длиной 1 млрд. заявок, на два порядка превышающей суммарную длину (9,01 млн заявок) выполненных на GPSS независимых прогонов, для w получена оценка 114,2, имеющая погрешность 4,4 % (почти на два порядка худшую). Таким образом, приведенный ко времени выигрыш в точности оценки составил около 100(100)2 = 106 раз.

7. Метод ARAND: решение проблемы искажения РТХ в ИМ. В описанном выше эксперименте на GPSS транзакты генерировались блоком GENERATE (Pareto(1,1,1.1)). Естественно, при этом моменты распределения Парето F (t) = 1 — (K /t) существенно смещались [2]. Но это не сказалось на точности расчета w, поскольку в СМО Pa/M/1 величина w не зависит от моментов интервала поступления заявок, а определяется его преобразованием Лапласа. Точные расчеты (опускаемые здесь из-за недостатка места) показывают, что дискретизация 6-разрядным ГСЧ распределения времени между приходами заявок не приводит к изменению величины w в моделируемой системе Pa/M/1 в пределах первых пяти-шести значащих цифр.

Что касается оценок средней длины очереди L и коэффициента загрузки р, то они оказываются сильно смещенными. Поскольку ГСЧ Uniform выдает в GPSS 6-разрядные числа, то реализуемое функцией Pareto(1,1,1.1) время между приходами заявок имеет среднее 8,0297, а не 11 [6]. Соответственно, оценка для р сходится в рассматриваемой СМО не к точному значению 4,4/11 = 0,4, а к величине 4,4/8,0297 = 0,548 (в описанном выше эксперименте на GPSS оценка для р получилась равной 0,546). При моделировании этой системы на AnyLogic, где ГСЧ реализует 15-разрядные числа, оценка для р сходится к 4,4/10,549 » 0,417. Это значение смещено относительно 0,4 на 4,3 %.

В [7] для решения проблемы искажений РТХ в ИМ предложен метод генерации с.в. с РТХ, названный каскадным. Приведем наиболее простую и точную его основу, не указанную в [7], — метод ARAND (Accurate RAND) генерации стандартных случайных чисел, представленный на рис. 4.

Метод ARAND преобразует обычные л-разрядные стандартные случайные числа (ССЧ)

Рис. 3. Усредненный ПП оценки для w в СМО Pa/M/1 при X = 1/11, р = 0,4

Рис. 4. Схема метода ARAND

z' в такие ССЧ z, которые, сколь бы малы они ни были, имеют n значащих цифр. Благодаря этому, с.в. с РТХ, реализуемые путем обратного преобразования хвоста распределения, тоже имеют n точных значащих цифр, и проблема искажения РТХ устраняется (подтверждающие это расчеты приведены в [7]). Например, обращая хвост распределения

Парето F (t) = 1 — F(t) = (K /1)a , получаем для генерации с.в. x, распределенной по Парето, формулу x = Kz~1/а , и, если z генерируется методом ARAND, то моменты с.в. x не смещаются. На рис. 5 приведена процедура, реализующая для GPSS World метод ARAND на языке Plus. Здесь Arg — это номер ГСЧ Uniform.

Первые 10 чисел на выходе ГСЧ Uniform(1,0,1) таковы: 0.842366, 0.777717, 0.880991, 0.260493, 0.463553, 0.083898, 0.022383, 0.948592, 0.344245, 0.929906. При обращениях к процедуре ARAND(1) они преобразуются в числа 0.842366, 0.777717, 0.880991, 0.260493, 0.463553, 0.00948592, 0.344245, 0.929906. Из-за потери части чисел исходной последовательности z' длина периода у ГСЧ ARAND в среднем в 1/0,9 раз меньше, чем у исходного ГСЧ. Период ГСЧ ARAND(1) содержит 1843200000 чисел, наименьшее из них равно 6.54026e-11.

Обращение к процедуре ARAND( ) осуществляется так же, как и к другим PLUS-процедурам [1]. Например, для генерации потока транзактов с временем между ними, распределенным по Парето с параметрами K = 1, а = 1,1, можно записать блок GENERATE следующим образом:

GENERATE (1#ARAND(1)A(-1/1.1))

Описанный в п. 6 эксперимент на GPSS в случае такого использования ГСЧ ARAND дает надежный степенной прогноз стационарного р, равный 0,397, ошибка которого по отношению к точному значению р = 0,4 составляет 0,8 %. Таким образом, с помощью процедуры ARAND (6-разрядной) получена оценка почти на порядок более точная, чем предельно достижимая оценка 0,417, обеспечиваемая (например, в AnyLogic) стандартным 15-разрядным ГСЧ.

8. Выводы. Модельное время в GPSS World иногда течет вспять, т.е. событиям-следствиям, имитируемым после наступления событий-причин, приписываются меньшие значения времени, чем событиям-причинам. Соответствующие ошибки оценок ВВХ растут во времени с большим ускорением и становятся неприемлемыми при длине прогона модели СМО более 100 млн заявок. При этом генераторы класса Uniform, с точки зрения равномерности распределения и статистической неза-

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

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

Проблемы искажения РТХ, возникающие из-за ограниченной разрядности стандартных ГСЧ, полностью устраняются при использовании предложенного в статье метода ARAND и процедуры, реализующей его в GPSS World.

Библиографический список

1. Руководство пользователя по GPSS World : пер. с англ. / Под ред. К. В. Кудашова. — Казань : Мастер Лайн, 2002. — 384 с.

2. Задорожный, В. Н. Проблемы генерации случайных величин с фрактальными распределениями / В. Н. Задорожный, О. И. Кутузов // Омский научный вестник. Сер. Приборы, машины и технологии. - 2012. - № 3 (113) - С. 20 24.

3. Задорожный, В. Н. О качестве программных генераторов случайных чисел / В. Н. Задорожный // Омский научный вестник. Сер. Приборы, машины и технологии. - 2009. -№ 2 (80). - С. 199-205.

4. Карпов, Ю. Г. Имитационное моделирование систем. Введение в моделирование с AnyLogic 5 / Ю. Г. Карпов. -СПб. : БХВ-Петербург, 2005. - 400 с.

5. Zadorozhnyi, V. N. Simulation modeling of fractal queues / V. N. Zadorozhnyi // Dynamics of Systems, Mechanisms and Machines. Dynamics 2014. - Proceedings, 2015. PROCEEDINGS RUSSIA, OMSK, NOVEMBER 11-13, 2014 Р. 1-4

6. Задорожный, В. Н. Проблемы и техника моделирования фрактальных очередей / В. Н. Задорожный, О. И. Кутузов // Имитационное моделирование. Теория и практика : материалы 6-й Всерос. конф. В 2 т. Т. 1. - Казань : ФЭН, Академия наук РТ, 2013. - С. 143-148.

7. Задорожный, В. Н. Каскадный метод реализации распределений с тяжелыми хвостами / В. Н. Задорожный. - Омский научный вестник. Сер. Приборы, машины и технологии. -2015. - № 2 (140). - С. 222-226.

ЗАДОРОЖНЫЙ Владимир Николаевич, доктор технических наук, доцент (Россия), профессор кафедры автоматизированных систем обработки информации и управления. Адрес для переписки: [email protected]

Статья поступила в редакцию 15.09.2015 г. © В. Н. Задорожный