CC BY
95
\ ИНФОРМАЦИОННЫЕ КАНАЛЫ И СРЕДЫ
УДК 519.872
РАСЧЕТ СИСТЕМ ОБСЛУЖИВАНИЯ С ГРУППОВЫМ ПОСТУПЛЕНИЕМ ЗАЯВОК
Ю. И. Рыжиков,
доктор техн. наук, профессор
Военно-космическая академия им. А. Ф. Можайского
Предлагается способ расчета распределения времени ожидания заявок в одно- и n-канальной системе при простейшем потоке пачек заявок случайного объема с ограниченным размахом, а также распределения числа заявок в системе. Точность расчетов иллюстрируется сопоставлением с результатами имитационного моделирования.
A method is proposed to compute the number-in-the-system and waiting time distributions for a one- and multichannel queuing systems with a Poissonian flow of the random (finite volume) demands. The results are compared with imitation ones.
Введение
Работа многих систем обслуживания связана с обслуживанием неординарного потока заявок, поступающих пачками случайного, вообще говоря, объема. Такие ситуации порождаются прежде всего расщеплением первичных заявок, например при передаче данных в сетях с коммутацией пакетов, при заявках на комплекты изделий. Другими примерами могут служить прибытие групп пассажиров (семейных, экскурсионных, командированных) в пункт отправления, медицинская помощь жертвам катастроф и террористических актов, прорыв группой самолетов или ракет зоны ПВО и т. п.
Групповое поступление заявок существенно ухудшает показатели системы в сравнении с ординарным потоком той же средней интенсивности. Для оценки работы системы и последующего принятия организационных мер (управления потоком, увеличения быстродействия обслуживающих устройств или их числа) необходимо уметь рассчитывать следующие распределения:
— времени ожидания начала обслуживания пачки;
— дополнительных задержек внутри пачки и средней задержки заявок пачки;
— задержки пачки в целом;
— числа заявок, находящихся в системе.
Поставленная задача уже довольно давно обсуждается в литературе [8-12, 14]. Однако:
— все эти результаты относятся только к одноканальным системам;
— как правило, они имеют частный характер [14], а в ряде случаев [11] ошибочны;
— верификация результатов и данные о численной реализации отсутствуют.
В данной статье последовательно рассматриваются подходы к расчету систем обслуживания групповых заявок, решающие упомянутые проблемы. Объем пачек предполагается ограниченным, а входящий поток пачек — простейшим.
Имитационное моделирование обслуживания групповых заявок
Имитационное моделирование (ИМ) в рассматриваемых ситуациях в связи с отсутствием альтернативных (и даже просто апробированных) методов расчета является практически единственным способом получения эталонных результатов. В инструментальных системах ИМ, например в GPSS World [1, 4, 7], предусмотрены ситуации с расщеплением и последующей сборкой заявок, означающей окончание обслуживания пачки. Однако ответа на все поставленные вопросы любая система со встроенным интерпретатором не дает. В GPSS World нельзя, к примеру, непосредственно получить моменты распределения времени ожидания пачки порядка выше второго, распределение времени пребывания в системе в зависимости от номера заявки в пачке, распределение числа заявок в системе. Встроенный в упомянутую систему язык PLUS все же является усеченным подмножеством универсальных алгоритмических языков, ориентированных на численные приложения, и по
возможностям и удобству заметно уступает Фортрану (см. раздел по оценке GPSS World [4]). Методика и приемы построения имитационных моделей на Фортране обсуждались в работах [4, 5]. Для решения вышеперечисленных задач имитационная модель должна иметь следующую специфику.
• По прибытии заявки случайным образом (в соответствии с заданным распределением) формируется объем пачки. Соответственно определяются количества d — принимаемых в каналы заявок и l — направляемых в очередь. Для первых свободные каналы определяются в цикле просмотра моментов освобождения с выходом из него после выявления d каналов. Свободные каналы занимаются, для них формируются моменты освобождения; идет подсчет количества принятых заявок.
• Для заявок, направляемых в очередь, определяется возможность их приема (очередь ограничена). Для всех принимаемых запоминается момент их прибытия. Кроме того, для головной заявки пачки в ее паспорте фиксируется объем пачки, а для последующих — ее номер в пачке. В случае приема на обслуживание хотя бы одной заявки пачки в паспортах остальных запоминается время начала обслуживания пачки.
• При завершении обслуживания и наличии заявок в очереди прежде всего выясняется статус головной заявки очереди. Если она — первая в своей пачке, то для остальных заявок пачки фиксируется момент начала ее обслуживания. Далее накапливаются степени истекшего времени ожидания. Если заявка — не первая, копятся степени дополнительной задержки по отношению к первой заявке пачки. Накопление идет раздельно по номерам заявок в пачке.
• Для построения распределения числа заявок в очереди фиксируются моменты его изменения, связанные с постановкой в очередь или началом обслуживания первой заявки пачки (последнее необходимо для расчета распределения ожидания начала обслуживания). Соответственно, имеется счетчик текущего числа пачек, и при его изменении к соответствующим накопительным ячейкам добавляется интервал неизменности.
• Время ожидания начала обслуживания пачки определяется как разность между временем выбора на обслуживание первой заявки пачки и моментом ее прибытия в систему. Суммирование интервалов и количества реализаций идет в раздельных (по номерам заявок в пачке) счетчиках, причем количество реализаций увеличивается и при немедленном приеме заявки на обслуживание.
Верификация модели проводилась сопоставлением результатов ее работы с полученным аналитическим расчетом системы MX /G/1 (см. далее).
Распределение числа заявок, прибывших за время обслуживания
Рассмотрим пуассоновский поток пачек требований, каждая из которых имеет одно и то же рас-
пределение числа заявок. Тогда распределение
{&п*} суммарного числа заявок в п пачках будет
п-кратной сверткой распределения {//*}={} .
Для дальнейшего нам необходимы вероятности прибытия ровно і заявок за случайный интервал времени между смежными обслуживаниями, имеющий распределение В(і):
h =J є'1' t ^ffdB(t)=
J “ П I
n=0
= Е г* / * (1)
Интегралы следует вычислить предварительно согласно рекомендациям разд. 3.6.3 работы [3], а свертки (для конечного размаха) последовательно получать численно и выполнять подсум-мирование раздельно для каждого I. Нужно иметь в виду происходящее на каждом шаге свертки удлинение массива вероятностей {/¿и*}.
В табл. 1 приводится сопоставление расчета по вышеописанной схеме и моделирования (1 млн испытаний) числа заявок обобщенного пуассонов-ского потока.
Объем пачки предполагался равновероятным в диапазоне 1^6, интервал времени — равномерно распределенным на интервале [0, 10].
Расчет одноканальной системы
Относительно стандартной системы М/в/1 известно (см., например, работу [3]), что стационарные вероятности наличия в системе ровно Ь заявок вычисляются как
Ро = 1 -И; (2)
( ь-1 Л
Pk =
Pk-1 - PoVk-1 -t PiVk-j j=1
/%, k = 1,2,...,
где b1 — среднее время обслуживания заявки, а
q=J e-1 MdBW
0
есть вероятность прибытия ровно заявок за случайное время обслуживания с распределением В(Ь). Преобразование Лапласа—Стилтьеса (ПЛС) распределения времени ожидания получается согласно известной формуле Полячека—Хинчина (ФПХ)
т(в) =---г-Ро--------------------, (3)
1 ~ [1 -Р(в>]
где Р(в) есть ПЛС от плотности В^).
■ Таблица 1. Обобщенный пуассоновский поток
п Расчет Модель п Расчет Модель
0 .43143е+00 .43237е+00 20 .16770е-02 .17110е-02
1 .49974е-01 .49654е-01 21 .12536е-02 .12640е-02
2 .54396е-01 .53787е-01 22 .91291е-03 .88300е-03
3 .59148е-01 .59052е-01 23 .64986е-03 .67300е-03
4 .64250е-01 .64103е-01 24 .45644е-03 .50100е-03
5 .69724е-01 .69234е-01 25 .32116е-03 .32400е-03
6 .75594е-01 .75437е-01 26 .22857е-03 .24000е-03
7 .31910е-01 .32270е-01 27 .15886е-03 .18600е-03
8 .29800е-01 .30080е-01 28 .10808е-03 .12900е-03
9 .27174е-01 .27388е-01 29 .72346е-04 .91000е-04
10 .23976е-01 .24203е-01 30 .48001е-04 .55000е-04
11 .20146е-01 .20224е-01 31 .31759е-04 .46000е-04
12 .15622е-01 .15626е-01 32 .20907е-04 .30000е-04
13 .10335е-01 .10292е-01 33 .13506е-04 .19000е-04
14 .86349е-02 .86930е-02 34 .86336е-05 .12000е-04
15 .70107е-02 .68830е-02 35 .55360е-05 .80000е-05
16 .55038е-02 .53480е-02 36 .36252е-05 .60000е-05
17 .41612е-02 .40850е-02 37 .24694е-05 .30000е-05
18 .30353е-02 .29600е-02 38 .17744е-05 .10000е-05
19 .21854е-02 .21270е-02 39 .13579е-05 .00000е+00
таблицы у(в) в окрестности нуля с последующей сменой знака у нечетных производных или численной сверткой в моментах на основе символического равенства ск = (а + Ь)к с переносом показателей степеней в индексы. Обе эти технологии при разумном выборе шага построения и длины таблицы дают практически совпадающие результаты и могут применяться при решении других задач. Тогда моменты распределения времени ожидания начала обслуживания пачки могут быть найдены по формулам (4) с заменой {Ь} на {^}. Вероятность свободного состояния системы здесь вычисляется как
Ро = 1 -АА, (6)
где ^ — средний объем пачки. Последующие вероятности распределения числа пачек вычисляются по формулам, аналогичным (2), с вычислением {дг} для определяемого (5) распределения трудоемкости пачки.
Распределение суммарного количества заявок в системе с групповым потоком следует рассчитывать согласно системе (2) с заменой первой формулы на (6) и вероятностей {д;} на вычисляемые согласно (1)вероятности{Н}.
Приведем сравнительные результаты моделирования (1 млн пачек равновероятного объема от 1 до 6) и расчета системы с равномерной на [0, 10]
Избавившись от знаменателя в правой части формулы (3), разложим входящие в нее преобразования Лапласа по степеням в и приравняем коэффициенты при одинаковых степенях в правой и левой частях. Из этих равенств следуют формула для первого момента распределения длительности ожидания пачки и рекуррентные формулы для высших моментов:
ХЪ2 »і =----------—
и —---------»
1 2(1 -ХЪ1)
юъ =7~ТЪТ^ -иъ +1 л1Ък+1-1ю1, ъ = 2,3■", (4) 1 -Щ Ро]!(Ъ + 1 -у)! ] ]
при этом ю0 = 1.
В одноканальной системе заявки пачки обслуживаются строго последовательно, а ПЛС времени обработки каждой выражается через ПЛС трудоемкости отдельной заявки Р(в) формулой
м
тс«) = £ т [рооГ=т*)), (5)
т=1
где F(•) — производящая функция объема пачки. Моменты {] этой трудоемкости можно получить многократным численным дифференцированием
■ Таблица 2. Распределение числа пачек в очереди по Ы/О/1
] Расчет Модель ] Расчет Модель
0 .40341е+00 .40224е+00 10 .85699е-02 .88636е-02
1 .16340е+00 .16312е+00 11 .61429е-02 .63665е-02
2 .12174е+00 .12099е+00 12 .44032е-02 .44870е-02
3 .88167е-01 .88313е-01 13 .31563е-02 .31766е-02
4 .63288е-01 .63167е-01 14 .22624е-02 .21522е-02
5 .45342е-01 .45830е-01 15 .16217е-02 .16070е-02
6 .32483е-01 .32488е-01 16 .11625е-02 .11304е-02
7 .23276е-01 .23682е-01 17 .83326е-03 .78494е-03
8 .16681е-01 .17140е-01 18 .59728е-03 .60351е-03
9 .11953е-01 .12123е-01 19 - .50317е-03
длительностью обслуживания заявки для коэффициента загрузки 0,8 (табл. 2).
Кроме того, были получены моменты длительности ожидания для пачки (табл. 3).
Здесь ФПХ подразумевает расчет высших моментов по рекуррентным формулам (4), а MFACT — через факториальные моменты {т[к]} длины очереди пачек согласно формуле Брюмелля [3]
к = т[к]/ Ак.
■ Таблица 4. Распределение числа заявок в системе Мх/С/1
] Расчет Модель ] Расчет Модель
0 .20025е+00 .19985е+00 21 .11847е-01 .11988е-01
1 .49231е-01 .49081е-01 22 .10809е-01 .10988е-01
2 .52554е-01 .52130е-01 23 .98625е-02 .99888е-02
3 .54630е-01 .54487е-01 24 .89987е-02 .91151е-02
4 .55031е-01 .54938е-01 25 .82106е-02 .83720е-02
5 .53278е-01 .52955е-01 26 .74914е-02 .76478е-02
6 .48852е-01 .48722е-01 27 .68353е-02 .69639е-02
7 .41197е-01 .41141е-01 28 .62366е-02 .63042е-02
8 .38521е-01 .38470е-01 29 .56904е-02 .57371е-02
9 .35595е-01 .35358е-01 30 .51920е-02 .52330е-02
10 .32594е-01 .32395е-01 31 .47372е-02 .47071е-02
11 .29685е-01 .29406е-01 32 .43223е-02 .43543е-02
12 .27006е-01 .26686е-01 33 .39438е-02 .39469е-02
13 .24627е-01 .24475е-01 34 .35984е-02 .36217е-02
14 .22507е-01 .22442е-01 35 .32832е-02 .33337е-02
15 .20541е-01 .20569е-01 36 .29956е-02 .30517е-02
16 .18737е-01 .18869е-01 37 .27333е-02 .28815е-02
17 .17093е-01 .17107е-01 38 .24939е-02 .26155е-02
18 .15595е-01 .15729е-01 39 .22755е-02 .24088е-02
19 .14231е-01 .14347е-01 40 .20762е-02 .22400е-02
20 .12985е-01 .13017е-01 41 .18944е-02 .20520е-02
■ Таблица 3. Моменты распределения ожидания для пачки
Способ расчета Порядок момента
1 2 3
По ФПХ .46594е+02 .51851е+04 .86156е+06
Через МЕАСТ .46594е+02 .51851е+04 .86156е+06
В модели .46938е+02 .52437е+04 .86378е+06
В модели через МЕАСТ .46943е+02 .52475е+04 .86595е+06
Распределение числа заявок в системе (табл. 4), как отмечалось выше, можно получить по алгоритму для стандартной системы М/О/1 после замены {д-} на их аналоги для обобщенного пуассо-новского потока из табл. 1.
Обращает на себя внимание вызванное группировкой заявок в пачки сильнейшее затягивание «хвостов» распределения.
Поскольку дополнительная задержка і-й заявки пачки имеет ПЛС фі(в) = Рі-1(в), для пачки в целом ПЛС задержки
M
эт-1
(s).
Фе («) = Е /тР"
т=1
Для средней задержки произвольной заявки пачки имеем
__ M f m
Ф(«) =t — -1(s) =7-
m=1 mu 1
M
-Ґ1
[1 -Pm (s)].
--P(s) m=1m
Распределение полного времени пребывания заявки в системе получается сверткой распределений ожидания пачки, дополнительной задержки в пачке и чистой длительности обслуживания.
Многоканальная система
Прежде всего отметим, что наиболее удобным и универсальным методом расчета многоканальных систем с произвольным распределением обслуживания является аппроксимация последнего гипер-экспоненциальным Н2. Такая аппроксимация позволяет сохранить три момента исходного распределения, что можно считать необходимым и достаточным.
Ниже обсуждается модель с ординарным потоком заявок, на основе которой далее предлагается метод для «групповой» задачи.
Диаграммы переходов и расчет М/Н2/п
Пусть вероятность того, что длительность обслуживания превышает ' (она же ДФР — дополнительная функция распределения):
B(') = £ yie~^i'. i=1
Тогда задачу можно рассматривать как процесс обслуживания заявок двух типов, причем тип заявки назначается с вероятностями {yi} в момент выбора ее на обслуживание. Характеризуя состояние системы полным числом заявок в ней (номер яруса) и распределением обслуживаемых заявок по типам, получаем диаграммы переходов по прибытии заявок и завершению обслуживания — для трехканальной системы (рис. 1 и 2 соответственно).
На рис. 2 при j >n поток обслуживания заявок i-го типа равен qi ці, где qi — содержимое i-й позиции кода микросостояния. При наличии очереди завершение обслуживания в зависимости от типа
выбранной из очереди заявки с вероятностями {у} приводит в одно из двух микросостояний вышележащего яруса.
На основе диаграммы переходов могут быть получены матрицы интенсивностей переходов:
Лу — с у-го яруса на (у + 1)-й по прибытии заявки, Bj — с у-го на (у - 1)-й по завершению обслуживания,
Бу — ухода из микросостояний у-го яруса (диагональные).
В данной модели нет переходов в пределах одного яруса, так что фигурирующие в общем случае
1У1
1 0 1
1У1
3 о
3 о
О 3
О 3
Рис. 1. Переходы по прибытии заявки в системе М/Н2/3
О 0 0
4 3 0
2 1
1 2
0 3
0 3
Рис. 2. Переходы по обслуживанию заявки в системе М/Н2/3
1
1
[3] матрицы {Су} здесь отсутствуют. Микросостояния яруса-источника соответствуют строкам, а яруса-приемника - столбцам матриц. Любой ненулевой элемент матрицы равен метке на представляющей переход стрелке. Элементы матриц {Ву} для четвертого и последующих ярусов получаются умножением исходной («корневой») интенсивности на у, соответствующий состоянию-приемнику. Поскольку для у > п матрицы Ау = XI, Бу = = Бп, а для у > п + 1 Ву = Вп+!, запоминать нужно сравнительно небольшое число матриц.
С помощью этих матриц можно записать векторно-матричные уравнения баланса переходов между ярусами. Решение упомянутых уравнений можно получить методом матрично-геометрической прогрессии (см. [3]) или восходящим к Така-хаси и Таками итерационным методом [2, 13]. Последний далее развивается нами применительно к потоку групповых заявок.
Распределение количества пачек и моменты распределения ожидания
Распределение количества находящихся в системе Мх/Н2/п пачек заявок заманчиво рассчитывать аналогично одноканальному случаю, т. е. применяя Н2-аппроксимацию к распределению трудоемкости пачки. Далее моменты распределения ожидания пачки могут быть найдены по распределению числа пачек в очереди. Имитационные эксперименты подтвердили хорошее согласие полученных результатов с непосредственным определением ожидания головной заявки пачки.
Этот подход неявно предполагает, что все заявки каждой пачки последовательно обслуживаются в одном и том же канале. Тем самым исключается разделение пачки между несколькими каналами в недогруженной системе и снижается реальная производительность последней. Можно предполагать, что возникающая погрешность будет возрастать при увеличении числа каналов и уменьшении номинального коэффициента загрузки системы:
р=Х/Ь1/п. (7)
Упоминавшиеся ожидаемые тенденции иллюстрируются сопоставлением результатов расчета (Р) и имитации (И) в табл. 5.
Таким образом, описанный подход вполне может найти практическое применение, но в ограниченном диапазоне условий.
Дополнительные задержки в пачках
Эти задержки могут быть рассчитаны аналогично одноканальному случаю с заменой распределения времени на распределение интервалов между выбираемыми из очереди заявками. ДФР этого распределения можно [6, с. 253] аппроксимировать формулой
Вп (О =
В*(0
п-1.
В(0, (8)
которая имеет прозрачный вероятностный смысл: для одного из каналов (только что завершившего обслуживание) берется полное распределение длительности обслуживания, а для прочих — помеченное звездочкой остаточное (случайная модификация полного).
Возможность применения формулы (8) была проверена на имитационной модели, предварительно откалиброванной по задаче с показательным распределением длительности обслуживания (в этом случае остаточное распределение совпадает с исходным и Вп(¿) = е~щ{:). Рабочее моделирование выполнялось для распределений длительности обслуживания Эрланга 3-го порядка (Е3), коэффициент вариации V = 0,577, и гиперэкспо-ненциального Н2, V = 1,36, при одинаковой средней длительности Ь1 = 5. Результаты расчета моментов Вп (£) на основании формулы Саати (С) и посредством модели (М) представлены в табл. 6.
Анализ таблицы с учетом общеизвестного факта быстрого роста относительных погрешностей моделирования при увеличении порядка вычисляемых моментов указывает на допустимость практического использования формулы (8).
Еще одним вариантом описания потока обслу-живаний в полностью загруженной системе является суммирование потоков от п каналов по методике, описанной в работе [3, разд. 3.6].
Наконец, для решения задачи можно воспользоваться распределением Вейбулла. Здесь ДФР имеет вид В^) = ехр(-^/Т), а вероятность того, что обслуживание не завершится ни в одном из п каналов:
Вп (£) = ехр(-п^ /Г), (9)
т. е. описывается тем же распределением с заменой Т на Т/п. Подставляя это значение в формулы для моментов распределения Вейбулла, убеждаем-
■ Таблица 5. Моменты распределения ожидания пачек
Моменты © II 50 © =
п = 3 п = 5 п=3 п=5
И Р И Р И Р И Р
Ш1 .145е+2 .128е+2 .787е+1 .667е+1 .477е+1 .362е+1 .244е+1 .149е+1
Ш2 .550е+3 .483е+3 .170е+3 .152е+3 .838е+2 .674е+2 .265е+2 .173е+2
^3 .307е+5 .268е+5 .531е+4 .510е+4 .217е+4 .180е+4 .428е+3 .284е+3
■ Таблица 6. Моменты интервалов между обслуживаниями в п-канальной системе
Тип В(Ь) Порядок момента Число каналов п
2 3 4
М С М С М С
Е3 1 2.50 2.50 1.62 1.67 1.22 1.25
2 9.48 9.46 4.28 4.46 2.50 2.60
3 4.64е1 4.57е1 1.47е1 1.57е1 6.83 7.21
4 2.75е2 2.64е2 6.08е1 6.66е1 2.30е1 2.45е1
5 1.90е3 1.77е3 2.90е2 3.28е2 9.13е1 9.72е1
6 1.49е4 1.36е4 1.55е3 1.83е3 4.20е2 4.36е2
Н2 1 2.62 2.45 1.68 1.59 1.27 1.16
2 1.70е1 1.58е1 6.48 6.57 3.62 3.48
3 2.05е2 1.75е2 4.36е1 4.54е1 1.75е1 1.74е1
4 3.95е3 2.94е3 4.41е2 4.55е2 1.25е2 1.25е2
5 1.09е5 6.94е4 5.88е3 6.18е3 1.17е3 1.19е3
6 3.80е6 2.13е6 9.37е4 1.09е5 1.30е4 1.43е4
ся, что его моменты порядка т получаются делением исходных на пт/к. Итак, здесь надлежит:
— найти параметры аппроксимации распределением Вейбулла остаточного распределения длительности обслуживания по двум моментам;
— вышеупомянутым пересчетом (с заменой п на п -1) получить три момента распределения (9);
— аппроксимировать это распределение гипер-экспоненциальным;
— найти такую же аппроксимацию для исходного распределения;
— выполнить суммирование двух потоков с найденными распределениями интервалов между заявками.
Распределение числа заявок в системе
Мх/Н2/п
Запишем условия баланса для векторов-строк {уу} вероятностей микросостояний системы с учетом прибытия заявок пачками случайного объема не свыше М:
М т-1
Уі“ і = Х1 т У і-т П А і-т+І +У і+1В/+1- (10)
т=1 1=0
Здесь все матрицы {Аі} разделены на интенсивность X потока пачек. Каждое слагаемое с множителем іт соответствует пачке из т заявок, так что для попадания на у-й ярус исходный должен иметь номер і - т. Сомножители произведения с учетом вышеупомянутого деления имеют смысл вероятностей переходов по прибытии заявки между смежными ярусами, а все произведение — вероятности перехода между микросостояниями ярусов І - т и у. Предельный индекс суммирования
М = шт{у, М}, так что для нулевого яруса правая часть(10)сводится ку1В1.
Выразим векторы {уу} через векторы ^у} условных вероятностей микросостояний яруса и суммарные вероятности {р} наличия в системе ровно у заявок: уу = руЬу, и введем отношения вероятностей смежных ярусов
гу = Ру -1/ Ру, ху = Ру+1/ Ру. (11)
Теперь уравнения (10) можно переписать в виде
М (т-1 \( т-1 \
і-т
т=1
і-т+І
Таким образом:
* і = і'і + х$,
где
р' = х
М
( т-1 Л( т-1
1 У!*-,1' П і П Аі-т
т=1
1=1
+ хі * і+1В і+1.
(12)
(13)
О-1; (14)
1=0
р' = *і+> В і+1“ -!.
(15)
Верхними индексами (в скобках) указаны номера последовательных приближений, в каждом из которых делается прогонка по всем ярусам. Смысл перехода к векторам условных вероятностей в том, что только для них удается задать хорошие начальные приближения, причем лишь для больших у. Последнее обстоятельство определяет и направление прогонки — от больших значений индексов у к нулю. Второе преимущество перехода
к условным векторам вероятностей — возможность замыкания расчетной схемы на основе допущения
t(k) = t(k-i)
/max /max -1 '
Переход к условным вероятностям связан с введением для каждого яруса дополнительных неизвестных {Xj} и {Zj}. Соответственно нужны дополнительные уравнения для их определения. В качестве первого выберем баланс переходов между j-м и (/' +1)-м ярусами. Интенсивность переходов сверху должна учитывать М = min{M -1,/} ярусов выше j-го; она составит
Лi = Х
(
Pj + Pj-і (1 - /і)+Pj-2 а - а - /2)+• • •+Pj-м+1
i-1
1 -S /і
i=1
=Pj Х
i=1
1 + zj Si 1 -S /т П
m=1
Интенсивность переходов снизу есть Ру+1у)1Ву+11у. Приравнивая эти количества и разделив обе части равенства на Ру, получаем уравнение
Х + Хг;
i=1
i-1
M ( i
S 1 -S fm і П 2I -
А
m=1
m=1
= x t(k) B 1 xjl/+1R j+11/'
(16)
Второе дополнительное уравнение (гуру + +ХуРу)1у = 1 дает условие нормировки компонент 1у. Из него следует, что
Ху =(1 - ^ур'у 1у )/(у ). (17)
В этих формулах 1 у есть вектор-столбец из единиц с числом элементов, равным количеству микросостояний на у-м ярусе. Умножение на него справа соответствует подсчету суммы компонент левого сомножителя. Подставляя (17) в (16), убеждаемся, что
2; =-
в j+1lj/(( lj )-Х
M (
А( i-1
S 1 -S /т П
j-т
(18)
Выведенные соотношения для у=ушах, ушах-1, — , 1 применяются в следующей очередности:
1) получить ву и ру согласно (14) и (15);
2) найти 2у по формуле (18);
3) вычислить Ху согласно (17);
4) получить очередное приближение у по формуле (13).
Поскольку на нулевом ярусе имеется всего одно микросостояние, 10 = 1 и в пересчете не нуждается.
Начальные приближения
Для расчета начальных приближений воспользуемся установленной экспериментально стабилизацией векторов {1у} и отношений смежных вероятностей при у —— ^. В ряде работ по теории очередей для предельного Х предлагается формула
x = р2М +vB),
(19)
где р — коэффициент загрузки системы; ьА и ьВ — квадраты коэффициентов вариации распределений интервалов между прибывающими заявками и длительности обслуживания соответственно. В нашем случае коэффициент загрузки вычисляется согласно (7), а интервалы между заявками для второй и последующих заявок пачки равны нулю. Таким образом, для простейшего потока пачек интенсивности X и пачки объема т средний интервал между заявками
а1(т) = —1 + ^—1 • 0 = 1/(тХ), т X т
а второй момент а2(т) = 2/(тХ)2. Подставляя в выражение для квадрата коэффициента вариации = а2/ а^ -1 моменты интервалов между заявками, усредненные по распределению объема пачки, находим
( M
S /т /
т=1
А
т
-1.
(20)
Для определения предельного вектора 1 воспользуемся отмеченной выше стабилизацией при у > п матриц интенсивностей переходов и самих векторов {1у}. Тогда индексы при векторах, матрицах, {ху} и {гу} можно опустить, а все матрицы А, упоминаемые в уравнении (12), считать единичными и из рассмотрения исключить. Теперь уравнение (12) запишется в виде
tD = Х
( м
S /т
т=1
А
t + xtB,
(21)
где г = 1/х. Расписав это векторно-матричное уравнение покомпонентно и заменив одно из уравнений на условие нормировки компонент вектора 1, получаем систему линейных алгебраических уравнений
для начальных приближений к {1у}, у = п, ушах. Для
у = 1, п -1 компоненты начальных векторов {1у} приходится принять равновероятными.
Начальные значения {гу} для у > п следует считать как 1/х, где х определяется согласно (19) с подстановкой vA из (20). Начальные {гу} для первых ярусов диаграммы можно получить с помощью соображений, использованных при выводе (16). Проводя горизонтальные разрезы между ярусами у и у - 1, убеждаемся, что
Р(0) ==-[ Ру-1 + (1 - /1) Ру -2 + ••• + (1 - /1 - /м-1) Ру-М+1 ].
^у
Набор учитываемых вероятностей для каждого разреза ограничивается требованием неотрицательности индексов. Усредненные интенсивности обслуживания ру = у/Ь1. Вероятность р00) предполагается равной единице. Далее вычисляются начальные приближения к отношениям вероятностей
г(0) = Р(0}/Р<0), у = М.
Условие прекращения итераций
Итерации можно считать законченными, когда стабилизируются значения {гу}, т. е. будет вы-
/ {| z(k) - z(k 1) |} < е. Требова-
полнено условие max
ния по точности естественно задать на уровне практически значимых вероятностей: е = 10-6.
Расчет итоговых вероятностей состояний
Условие баланса заявок для я-канальной системы в данном случае имеет вид
я-1 _
X (я - j)Р/ = я-^/&1. (22)
/=0
Очевиден вероятностный смысл условия: ожидаемое число свободных каналов равно недогрузке системы. Входящие в (22) вероятности последовательно выражаются через p0: p1 = p0/z1, p2 = = p1/z2 = p0/(z1z2) и т. д. Следовательно:
Ро
я-1 / j
я+X(я -/у П z
/=1 / г=1
= я - Xfb1,
откуда и определяется Р0. Последующие вероят-ностиРу = Ру-1/гу, у = 1, 2,... .
Результаты расчета
Расчет по описанной методике выполнялся для трехканальной системы с равномерным на [0,10] распределением обслуживания и равновероятным (от 1 до 6) объемом пачек заявок. Интенсивность потока пачек выбиралась для получения коэффи-
циента загрузки р = 0,8. Опорное отношение вероятностей z^ = 1/x^ составило 1,111. Для jmax = 39 стабилизация {z/} (табл. 7) с точностью 10-9 потребовала 44 итераций.
Дальнейшие {z/} практически совпадают с z27. Как видно, они очень близки к опорному. То же можно сказать и о предельных вероятностях условных микросостояний.
В табл. 8 расчетное распределение числа заявок в системе сопоставлено с результатами ИМ (10 млн испытаний).
Результаты согласуются очень хорошо. Отметим также, что сумма полученных вероятностей, дополненная суммой образующих геометрическую прогрессию неучтенных членов, с точностью 10-7 совпала с единицей, причем условие нормировки полных вероятностей в расчетной схеме не использовалось. Это согласие является дополнительным подтверждением правильности расчетной схемы и ее программной реализации.
Имитационная модель допускала до 100 заявок в системе. Однако сходимость численного метода при увеличении jmax заметно ухудшается. Поэтому целесообразно выполнять расчет при умеренных jmax, а требуемое число дополнительных вероятностей вычислять по формуле Р; = P;--|/Z; .
J J •'max
Распределение пачек и их задержек
Выше отмечалась необходимость расчета распределения числа «непочатых» пачек в очереди. Такая задача является обратной по отношению к естественному переходу от распределения числа пачек к распределению числа заявок и как таковая обладает вычислительной неустойчивостью. Это подтвердили попытки решения соответствующего уравнения в дискретных свертках несколькими разными методами. В связи с этим ниже предлагается хорошо зарекомендовавший себя полу-эмпирический метод.
Анализ результатов имитации показал, что отношение смежных вероятностей искомого распределения близко к интуитивно ожидаемому: Z = Zf. Полные пачки в очереди отсутствуют, если заявок в очереди нет или их количество^ превышает среднего остатка пачки (скажем, k — поло-
■ Таблица 7. Расчетные отношения смежных вероятностей
j zj j z j zj j Z
0 .10000e+01 7 .11719e+01 14 .10929e+01 21 .10957e+01
1 .20589e+01 8 .10837e+01 15 .10940e+01 22 .10958e+01
2 .11822e+01 9 .10776e+01 16 .10952e+01 23 .10958e+01
3 .10695e+01 10 .10824e+01 17 .10958e+01 24 .10958e+01
4 .96669e+00 11 .10918e+01 18 .10958e+01 25 .10958e+01
5 .99908e+00 12 .10972e+01 19 .10956e+01 26 .10958e+01
6 .10468e+01 13 .10976e+01 20 .10956e+01 27 .10959e+01
■ Таблица 8. Распределения количества заявок в системе
у' Расчет Имитация у' Расчет Имитация
0 .13692е+00 .13727е+00 20 .14154е-01 .14266е-01
1 .66499е-01 .65292е-01 21 .12917е-01 .13041е-01
2 .56251е-01 .58327е-01 22 .11788е-01 .11928е-01
3 .52598е-01 .52512е-01 23 .10757е-01 .10876е-01
4 .54410е-01 .54403е-01 24 .98164е-02 .98637е-02
5 .54460е-01 .54402е-01 25 .89579е-02 .90459е-02
6 .52027е-01 .51411е-01 26 .81745е-02 .82060е-02
7 .44397е-01 .44084е-01 27 .74594е-02 .74814е-02
8 .40967е-01 .40911е-01 28 .68069е-02 .67880е-02
9 .38016е-01 .37984е-01 29 .62115е-02 .61825е-02
10 .35122е-01 .35059е-01 30 .56681е-02 .56444е-02
11 .32169е-01 .32082е-01 31 .51723е-02 .51133е-02
12 .29320е-01 .29277е-01 32 .47198е-02 .46397е-02
13 .26713е-01 .26761е-01 33 .43070е-02 .42310е-02
14 .24442е-01 .24510е-01 34 .39302е-02 .38639е-02
15 .22341е-01 .22394е-01 35 .35864е-02 .35175е-02
16 .20400е-01 .20592е-01 36 .32726е-02 .32098е-02
17 .18617е-01 .18672е-01 37 .29864е-0 .29245е-02
18 .16990е-01 .17102е-01 38 .27251е-02 .26566е-02
19 .15507е-01 .15581е-01 39 .24867е-02 .24467е-02
вины ее среднего значения). Тогда можно считать, что вероятность этого события
п+к
40 = X Рь,
¿=0
а остальные вероятности образуют убывающую геометрическую прогрессию с начальным членом
41 и знаменателем Z. Замыкающую схему недостающую вероятность находим из условия нормировки: 41 = (1 - 1^) (1 -40).
В табл. 9 проводится сопоставление результатов расчета с имитацией.
Имея распределение числа непочатых пачек в очереди, можно по формуле Брюмелля выразить
■ Таблица 9. Распределения количества пачек в очереди
у' Расчет Имитация у Расчет Имитация
0 .47316е+00 .45458е+00 15 .16276е-02 .14718е-02
1 .14443е+00 .15136е+00 16 .11814е-02 .10227е-02
2 .10484е+00 .11139е+00 17 .85751е-03 .71961е-03
3 .76095е-01 .80331е-01 18 .62243е-03 .51405е-03
4 .55234е-01 .57622е-01 19 .45179е-03 .37014е-03
5 .40092е-01 .41471е-01 20 .32793е-03 .25248е-03
6 .29101е-01 .29731е-01 21 .23803е-03 .16101е-03
7 .21123е-01 .20981е-01 22 .17278е-03 .12383е-03
8 .15332е-01 .14868е-01 23 .12541е-03 .98421е-04
9 .11129е-01 .10666е-01 24 .91029е-04 .68672е-04
10 .80779е-02 .76942е-02 25 .66074е-04 .50771е-04
11 .58634е-02 .55641е-02 26 .47960е-04 .33102е-04
12 .42559е-02 .39279е-02 27 .34812е-04 .19618е-04
13 .30892е-02 .28407е-02 28 .25268е-04 .13579е-04
14 .22423е-02 .20432е-02 29 .18341е-04 .68985е-05
■ Таблица 10. Моменты распределения ожидания пачек
Способ определения w1 w2 w3
Расчет .14013e+02 .54106e+03 .31337e+05
Имитация непосредственно .14079e+02 .51615e+03 .28040e+05
Имитация через MFACT .14067e+02 .51369e+03 .27659e+05
моменты распределения времени ожидания начала обслуживания пачки (табл. 10). Заключение Широкий круг практически важных систем обслуживания описывается математическими моделями с групповым потоком заявок. Предло- женные расчетные зависимости впервые дают возможность получать основные вероятностные характеристики многоканальных систем. Имитационные эксперименты подтвердили достаточную точность расчета практически значимых вероятностей и моментов временных распределений.
Литература /
1. Кудрявцев Е. М. GPSS World. Основы имитационного моделирования различных систем. М.: ДМК Пресс, 2004. 320 с. 2. Рыжиков Ю. И., Хомоненко А. Д. Итеративный метод расчета многоканальных систем с произвольным распределением времени обслуживания //Проблемы управления и теории информации. 1980. № 3. С.32-38. 3. Рыжиков Ю. И. Теория очередей и управление запасами. СПб.: Питер, 2001. 384 с. 4. Рыжиков Ю. И. Имитационное моделирование. Теория и технологии. СПб.: КОРОНА принт, 2004. 384 с. 5. Рыжиков Ю. И. Современный Фортран. СПб.: КОРОНА принт, 2004. 288 с. 6. Саати Т. Л. Элементы теории массового обслуживания и ее приложения: Пер. с англ. М.: Сов. радио, 1965. 510 с. 7. Томашевский В. Н., Жданова Е. Г. Имитационное моделирование в среде GPSS. М.: Бестселлер, 2003. 416 с. 8. Chaudry M. L., Gupta U. C., Agarwal M. Computational analysis of the distribution of demands number in the system МХ/G/1 — an alternative approach // INFOR. 1992. Vol. 30. N 1. P. 30-43. 9. Krakowski M. Arrival and departure processes in queues — Pollaczek — Khintchine formulas for bulk arrivals and bounded systems //RFAIRO. 1974. Vol. 1. P. 45-56. 10. Loris-Teghem J. A. Modeles d’attente M/G/1 et GI/ /M/1 а arrivees et service en groupe. Berlin: Springer, 1969. 53 p. 11. Ommeren J. C. W., van. Simple approximation for the batch-size MX/G/1 queue //OR. 1990. Vol. 38. N 4. P. 678-685. 12. Sivasamy R. A bulk service queue with accessible and non-accessible batches //Opsearh. 1990. Vol. 27. N 1. P.46-54. 13. Takahashi Y., Takami Y. A numerical method for the steady-state probabilities of a GI/G/c queuing system in a general class //J. of the Operat. res. soc. of Japan. 1976. Vol.19. N 2. P. 147-157. 14. Wirth K. D. A remark on relations between batch delays and customer delays // J. of Information Processing and Cybernetics. 1985. Vol. 21. N 1/2. P. 65-67.
/
