CC BY
147
УДК 519.872
ПОЛНЫЙ РАСЧЕТ СИСТЕМЫ ОБСЛУЖИВАНИЯ С РАСПРЕДЕЛЕНИЯМИ КОКСА
Ю. И. Рыжиков,
доктор техн. наук, профессор
Военно-космическая академия им. А. Ф. Можайского
Разработан алгоритм полного расчета стационарных характеристик системы обслуживания вида C2/C2/n: стационарного распределения числа заявок, их распределения перед прибытием очередной заявки, моментов распределения интервалов между обслуженными заявками. Обсуждаются частные случаи и обобщения. Рассмотрены особенности программной реализации алгоритма, приведены результаты тестирования. Алгоритм предлагается использовать при расчете сетей обслуживания методом потокоэквивалентной декомпозиции.
An algorithm is proposed to compute the characteristics of the system C2/C2/n. This algorithm computes the stationary distribution of the number of demands, its distribution before arriving, and the moments of interdeparting intervals. Particular cases and generalizations, programming realization and testing results are discussed. The algorithm is planned to be used in the flow-equivalent network analysis.
Введение
Известно значительное влияние на результаты расчета сложных систем и сетей обслуживания не только средних интервалов между заявками и длительностей обслуживания, но и их распределений. Универсальным методом получения требуемых числовых характеристик является имитационное моделирование, существенными недостатками которого являются сложная логика программ, большой расход машинного времени и низкая точность определения малых вероятностей. Системы автоматизации имитационного моделирования типа GPSS (включая новейшую модификацию последней - GPSS World) не снимают перечисленных проблем, а лишь упрощают программирование ценой замедления работы моделей (как показали проведенные автором эксперименты - в десятки раз). Таким образом, несмотря на увеличение быстродействия ЭВМ, численные методы расчета СМО сохраняют свою актуальность.
Опубликовано множество методов расчета различных классов немарковских систем обслуживания, наиболее универсальные из которых опираются на подбираемые по методу моментов фазовые аппроксимации исходных распределений показательным, эрланговским, гиперэкспоненциаль-ным, коксовым и общим распределением Ньютса. Для последнего, однако, отсутствуют алгоритмы расчета параметров аппроксимации. Эрланговские распределения из-за требования целочисленности фазового параметра не обеспечивают точного вы-
равнивания даже двух моментов. Поэтому наибольшую практическую ценность представляют гиперэкспоненциальное и коксово распределения. Однако для первого из них распределения Эрланга являются особыми случаями: система последовательных фаз не может быть сведена к тому же числу параллельных. Это гарантирует аварийное завершение программы подбора параметров аппроксимирующего распределения, вынуждает менять тип аппроксимации в зависимости от исходных данных и требует обширного набора программ, включающего все возможные комбинации. От данного недостатка свободно распределение Кокса, исследование свойств и систематическое использование которого в обеих позициях систем вида 01/0/п впервые описываются ниже.
Распределение Кокса
К фазовым распределениям относятся порождаемые системой подлежащих прохождению фаз обслуживания (прибытия заявки) с показательно распределенной длительностью пребывания в каждой из них. При фиксации номера фазы такие распределения приобретают марковское свойство, что и определяет целесообразность их использования. Для практических целей применение двухфазного распределения С2, обеспечивающего выравнивание трех начальных моментов, представляется вполне достаточным. Диаграмма этого распределения представлена на рис. 1. Параметры и |л,2 задают интенсивности потоков пе-
рехода, у - вероятность потребности во второй фазе, у - вероятность отсутствия таковой.
Получим алгоритм расчета параметров распределения С2 по его моментам {/} . Прежде всего заметим, что эти моменты с вероятностью у являются моментами свертки показательно распределенных задержек в фазах, а с дополнительной к ней суть моменты задержки в первой фазе. Таким образом, для расчета параметров С2 имеем систему уравнений
У_
М1
+ У
= /і;
У^+У Мі
у Мі М2 у с \
2 о 1 2
"Т + 2-----------+ _2'
Мі2 М1М2 м2
_ /2;
-6 6 „21 „12 6
У~3 + У _3 + 2---+ 3----2 +~3
Мі ч Мі Мі М2 М1 М2 М2 у
= /з-
Перейдем к обратным величинам |х; = М;1} и нормированным моментам / = / / І ! |. Тогда исходная система может быть переписана как
УХ1 + у(Х1 + Х2) = /1;
ух2 + у (х2 + х1х2 + 1 = /2;
ух3 + У (х3 + х2х2 + х1х2 + х3 | = /3.
С учетом у + у = 1 она преобразуется в систему х1 + ух2 = /1; х2 + ух2(х1 + х2) = х3 + ух2 (х2 + х1х2 + х% | = /3. (1)
Из первого уравнения системы исключаем
ух2 = /1 - х1 (2)
и подставляем в два оставшихся. Первое из них дает
_ /2 - х2 х2 ~ , х1.
/1 - х1
Конечная форма задачи - квадратное уравнение с решением
х1,2 _
/3 - /1/2 ±44/3/3 + 4/3 + /32 - 3/1 /2 - 6/1/2/3
2
(/ - /12 )
(3)
Решение для х1 берется с плюсом, у определяется из равенства (2). Заметим, что случай /2 = 1 соответствует одной экспоненте и должен выявляться на первом этапе алгоритма. Здесь есте-
■ Рис. 1. Двухфазное распределение Кокса
■ Таблица 1. Параметры коксовой аппроксимации
а у Ні Н2
0.2 -10.196 268 3.732
0.4 -3.682 .488 3.512
0.6 -1.584 .677 3.323
0.8 -.577 .845 3.155
1.0 .000 1.000 1.000
1.2 .745 2.853 1.147
1.4 .815 2.707 1.293
1.6 .880 2.555 1.445
1.8 .940 2.378 1.622
2.0 1.000 2.000 2.000
2.2 1.061-.011І 2.000+.354І 2.000-.354І
2.4 1.111-.027І 2.000+.485І 2.000-.485І
2.6 1.154-.044І 2.000+.577І 2.000-.577І
2.8 1.190-.062І 2.000+.649І 2.000-.649І
3.0 1.222-.079І 2.000+.707І 2.000-.707І
ственно принять у = 0 и ^ = 1/1\, а значение |Л,2 можно выбрать произвольно (например, равным ^,1).
Подставим в эти формулы моменты гамма-распределения с параметрами а и |1. Подкоренное выражение
^ а2 (а +1)(2-а)(а-1)2 °= 18/ .
Таким образом, при а > 2 аппроксимация должна иметь комплексные параметры. Случай а = 2 является граничным: он соответствует распределению Эрланга второго порядка с у = 1 и ^1 = |Л,2.
В табл. 1 приведены параметры ^-аппроксимации гамма-распределения при единичном первом моменте. Таблица подтверждает сделанный прогноз и снимает опасения относительно возможной отрицательности {|^}. Отметим парадоксальный, но ожидаемый результат: при а< 1 вероятность выхода процесса во вторую фазу отрицательна.
■ Таблица 2. Распределение числа заявок в системе Ы/О/1
] Параметр гамма-распределения
0.5 1.5 3.0 1е9
М01 МСЫ М01 МСЫ М01 МСЫ М01 МСЫ
0 .300е-0 .300е-0 .300е-0 .300е-0 .300е-0 .300е-0 .300е-0 .300е-0
1 .165е-0 .167е-0 .233е-0 .233е-0 .263е-0 .263е-0 .304е-0 .306е-0
2 .120е-0 .118е-0 .159е-0 .159е-0 .174е-0 .173е-0 .190е-0 .187е-0
3 .912е-1 .903е-1 .106е-0 .106е-0 .106е-0 .106е-0 .101е-0 .101е-0
4 .706е-1 .703е-1 .696е-1 .696е-1 .637е-1 .638е-1 .517е-1 .524е-1
5 .550е-1 .550е-1 .457е-1 .456е-1 .379е-1 .380е-1 .263е-1 .267е-1
6 .430е-1 .431е-1 .299е-1 .299е-1 .226е-1 .226е-1 .134е-1 .135е-1
7 .336е-1 .338е-1 .196е-1 .196е-1 .134е-1 .134е-1 .682е-2 .682е-2
8 .264е-1 .265е-1 .128е-1 .128е-1 .798е-2 .798е-2 .347е-2 .344е-2
9 .207е-1 .208е-1 .842е-2 .841е-2 .474е-2 .474е-2 .177е-2 .173е-2
10 .162е-1 .167е-1 .551е-2 .551е-2 .282е-2 .281е-2 .899е-3 .873е-3
11 .127е-1 .128е-1 .361е-2 .361е-2 .168е-2 .167е-2 .457е-3 .440е-3
12 .996е-2 .100е-1 .237е-2 .237е-2 .996е-3 .993е-3 .233е-3 .222е-3
13 .782е-2 .784е-2 .155е-2 .155е-2 .592е-3 .590е-3 .118е-3 .118е-3
14 .613е-2 .614е-2 .101е-2 .101е-2 .352е-3 .350е-3 .603е-4 .562е-4
15 .481е-2 .481е-2 .665е-3 .666е-3 .209е-3 .208е-3 .307е-4 .283е-4
16 .377е-2 .377е-2 .436е-3 .436е-3 .124е-3 .124е-3 .156е-4 .143е-4
17 .296е-2 .296е-2 .285е-3 .286е-3 .740е-4 .735е-4 .794е-5 .719е-5
18 .232е-2 .232е-2 .187е-3 .187е-3 .440е-4 .436е-4 .404е-5 .362е-5
Приведем результаты расчета распределения числа заявок в системе Ы/О/1 (табл. 2). Процедура MG1 использовала общеизвестные соотношения, реализующие метод вложенных цепей Маркова на базе рекуррентного расчета распределения {я} числа заявок, прибывающих за случайное время обслуживания:
для вырожденного и гамма-распределений [1]. Процедура MCN опиралась на С2-аппроксимацию и решала описанные ниже векторно-матричные уравнения баланса. Индекс ] указывает число заявок в системе. Параметр 1е9 фактически определяет вырожденное распределение, т. е. модель Ы/Б/1.
Гамма-распределение было выбрано как удобный эталон по двум причинам:
1) параметры его плотности
щ) = м(м^)а 1 е-м Г(а)
легко выражаются через среднее Ь1 и дисперсию Б согласно
а = Ь1 /Б, ц = а/Ь1;
2) вышеупомянутые для него вычисляются до удивления просто [1, с. 82]:
Я0 =
Я] = Я]-1
А, + |і X а + ] -1
^ + М ]
, ] = 1, 2, ... .
Хорошее согласие прямого расчета модели М/О/1 и конечных результатов ее коксовой аппроксимации в широком диапазоне условий, в том числе для комплексных и отрицательных параметров, убедительно свидетельствует о приемлемости обсуждаемой аппроксимации.
Достоинством С2-аппроксимации является естественное включение в нее как частных случаев М- и £2-моделей, а недостатками - трудность обобщения на большее число составляющих и усложнение диаграмм переходов.
Многофазное представление системы
Будем представлять состояние системы С2/С2/п полным числом заявок в системе, фазой прибытия очередной заявки и расстановкой проходящих обслуживание заявок по его фазам («ключами» соответствующих микросостояний). Например, ключ (21) означает, что две заявки проходят первую фазу обслуживания с интенсивностями ц1, а одна - вторую с интенсивностью ц2. Диаграмма состояний разделяется по фазам прибытия на две вертикальные ветви с интенсивностями Х1 и Х2 соответственно, а по числу заявок в системе - на горизонтальные ярусы. Совокупности ключей в левой и правой ветвях диаграммы дублируются. На рис. 2 показаны переходы по завершению фазы обслуживания для половины диаграммы (для другой половины они тождественны).
Обозначим через Б] множество всех возможных микросостояний системы, при которых на обслуживании находится ровно ] заявок, а через - количество элементов Б]. Нам потребуются следующие матрицы интенсивностей инфинитезималь-ных переходов:
А] ] ха]+1 ^ - в Б] + 1 (прибытие заявки),
в Б] (конец промежуточной фазы обслуживания или прибытия заявки),
В |>]х°]-1 ] - в Б]-1 (полное завершение обслуживания заявки),
Б] ] ха] ^ - ухода из состояний яруса ].
В квадратных скобках здесь и далее указывается размер матриц. Элемент (Ь, й) любой из этих матриц представляет интенсивность перехода из £-го состояния ]-го яруса в й-е состояние смежного (по переходам рассматриваемого типа) яруса. Заметим, что при ] > п все эти матрицы перестают зависеть от индекса. Ниже будут приведены примеры матриц для второго яруса.
Переходы по завершению фазы обслуживания не меняют фазы прибытия, что позволяет ограничиться изображением для них половины диаграммы (рис. 2). Здесь окончание первой фазы с интенсивностью у^1 переводит заявку во вторую фазу обслуживания, а с интенсивностью у^.1 - заверша-
00
Г\
^1у ^2
10 — ^у—► 01
|\ 1\
2^1у ^2 Ц1у 2^2
20 -2^1^^ 11 - М1у—► 02
ї\ 1\ 1\
3^1у ^2 2^1у 2^2 ^1у 3^2
30 — 3^у -► 21 — 2^у -► 12 — щу -► 03
1\ 1\ 1\
3^1у ^2 2^1у 2^2 ^1у 3^2
30 — 3^1у —^^21 — 2^1у —^^12 — ^1у
-03
\ \ \ \
Рис. 2. Интенсивности переходов по завершению обслуживания в фазе
ет обслуживание и поднимает изображающую точку на вышележащий ярус. Последнее при завершении второй фазы (интенсивность |Л,2) происходит всегда. Интенсивности переходов проставлены с учетом числа заявок, находящихся в каждой фазе исходного микросостояния.
Например, для полной диаграммы переходов по завершению обслуживания («удвоенный» рис. 2) имеем клеточно-диагональную матрицу интенсивностей переходов со второго на первый ярус:
Б2 =
2Міу 0 0 0
М2 М1у 0 0
0 2М2 0 0
0 0 2М1у 0
0 0 М2 М1у
0 0 0 2М2
Диаграмма переходов по фазам прибытия заявок в трехканальной системе оказывается весьма запутанной; поэтому мы ограничимся описанием логики ее построения. Если интервалы между заявками подчинены С2-распределению с параметра-
ми {и, Л15 Л2}, то все переходы из первой фазы с интенсивностями иХ1 приводят в аналогичное микросостояние второй фазы прибытия того же яруса, а с интенсивностями йХ1 - в первую фазу нижележащего яруса. Переходы по завершению второй (всегда последней) фазы имеют интенсивность Х2 и также приводят в первую фазу нижележащего яруса. Соответственно интенсивности «окончательных прибытий» описывает матрица
А2 =
иХ1 0 0 0 0 0 0 0
0 йХ1 0 0 0 0 0 0
0 0 их1 0 0 0 0 0
^2 0 0 0 0 0 0 0
0 ^2 0 0 0 0 0 0
0 0 ^2 0 0 0 0 0
Переходы в пределах яруса могут быть связаны с выходом из первой во вторую фазу обслуживания (в пределах каждой ветви) либо с прибытием заявки (переход из левой ветви в одноименное микросостояние правой):
Р =
0 2МіУ 0 иХ1 0 0
0 0 МіУ 0 иХ1 0
0 0 0 0 0 иХ1
0 0 0 0 2М 0
0 0 0 0 0 МіУ
0 0 0 0 0 0
Для диагональной матрицы Б2 суммарных интенсивностей ухода из микросостояний яруса мы приведем только ее диагональ:
diag Б2 = [Л1 + 2р,1, Х1 + р,1 + ц2,
^1 + 2^2, ^2 ^ 2^1, ^2 + ^1 ^^2, ^2 ^ 2^2 ].
Уравнения глобального баланса и их решение
Получим стационарные вероятности микросостояний, представленных на диаграммах предыдущего раздела.
Введем векторы-строки уу =<уу 1, уу 2, ..., у нахождения системы в состоянии (у, I), у = 0, 1, .... Теперь можно записать векторно-матричные уравнения баланса переходов между состояниями
Уо А =Уо Р +ііві;
Уі°і =уі-1аі-1 +УіРі +уі+1Б/+1, і =1 2, ••• • (4)
Система (4), дополненная условием нормировки, даже для моделей с ограниченной очередью
характеризуется чрезвычайно высокой размерностью, и стандартные методы решения систем линейных алгебраических уравнений применительно к ней оказываются малоэффективными. Для расчета разомкнутых систем обслуживания можно применить восходящий к Ивэнсу [5] метод матрично-геометрической прогрессии. Этот метод основан на представлении векторов вероятностей микросостояний для у > п формулой
У у =УпВ>~п, у = п п +1 ..., (5)
где Б - некоторая матрица (знаменатель прогрессии). Перепишем второе из уравнений (4) с учетом (5) для значения у, при котором произошла вышеупомянутая стабилизация матриц, так что их индексы можно опустить:
У у^ББ =У у _1 А + у у + у у _1Й2 В. (6)
Теперь ясно, что Б удовлетворяет матричному квадратному уравнению
Б2В + БЕ + А = 0, (7)
где Е = ¥-Б. Применим к его решению метод итераций. Будем искать поправку А к очередному приближению из уравнения (7), заменив в нем Б на Б + А. Пренебрегая членом, содержащим квадрат поправки, имеем
Б2В + АБВ + БАВ + БЕ + АЕ + А ~ 0
или
А(БВ + Е) + БАВ - -(Б2В + БЕ + А).
Последнее уравнение приводится к уравнению АО + БАВ = Н,
которое можно расписать как систему линейных алгебраических уравнений относительно компонент матричной поправки А и решить любым стандартным методом. В качестве начального приближения можно принять Б = р2/^А +1>В ^I, где р - коэффициент загрузки системы; иА и vB - коэффициенты вариации распределений интервалов между смежными заявками и длительности обслуживания соответственно; I - единичная матрица.
Далее усеченная система (4) для у = 0, п, в последнем из уравнений которой вектор уп +1 заменен на упБ, расписывается относительно компонент векторов {уу}, у = 0, п. Затем одно из уравнений заменяется на условие баланса заявок
п-1
X (п - /)У у 1у = п - Ь / а,
1-0
где 1у - вектор-столбец из а у единиц; а и Ь - средние интервал между заявками и длительность обслуживания. Наконец, находится решение системы для компонент начальных векторов.
Последующие векторы вероятностей микросостояний вычисляются согласно формуле (5).
Для решения системы (4) можно также применить менее быстродействующий, но более универсальный итерационный метод, впервые предложенный Такахаси и Таками [7] и описанный в работах [1, 2] в более общей форме, с детальной проработкой расчетных зависимостей и с частными вариантами.
Расчет выходящего потока
В сетях обслуживания потоки на входе узлов формируются из выходящих потоков других узлов. В общем случае выходящий поток сохраняет среднюю интенсивность (и средний интервал между заявками) входящего, но меняет вид распределения интервалов между заявками. Последнее удобно характеризовать коэффициентами немар-ковости
^ = d /d1 -i!, i = 2, 3, ...,
где {dj} - моменты распределения интервалов между заявками выходящего потока, который предполагается рекуррентным. Для простейшего потока упомянутые коэффициенты согласно данному определению равны нулю.
Ключом к проблеме преобразования потока в коксовом узле обслуживания является расчет марковской системы.
Марковские системы. В системе M/M/n, загруженной полностью, частное преобразование Лапласа-Стилтьеса (ПЛС) интервала до очередного обслуживания при наличии в системе j > n заявок
5j(s) = щ/(n^ + s), j = n, n +1, ... .
При j < n интенсивность обслуживания меняется с прибытием каждой очередной заявки. Будем интерпретировать s как параметр простейшего потока «катастроф». Тогда by (s) можно считать вероятностью отсутствия катастроф за интервал от ухода заявки, оставившей в системе ровно j заявок, до следующего завершения обслуживания. Вероятности появления событий каждого из рассматриваемых простейших потоков (катастроф, завершений обслуживания и прибытия новых заявок) пропорциональны их интенсивностям, т. е. s, j|j, и X соответственно. Теперь ясно, что
5; (s) =. ^ +. \ S/+i(s),
‘ j^ + A + s j^ + A + sJ
j = n -1, n - 2, ..., 0.
Здесь первое слагаемое соответствует завершению обслуживания, а второе - прибытию новой заявки (отчего число заявок в системе увеличилось на единицу) и отсутствию катастроф до следующего завершения обслуживания при новом начальном условии.
Результирующее распределение интервалов между заявками выходящего потока задается ПЛС
n-1
D(s)=Елі 5 і (s) + 1 і
n-1
5И (s),
где финальные вероятности {тсу} наличия в системе сразу после завершения обслуживания ровно у заявок при t подсчитываются согласно уравнению
i=1
(В)
Здесь {р} - стационарное распределение числа заявок в системе, а ^ равно /ц при <п и п^. - в противном случае. Знаменатель формулы (8) сводится к равенству
^ViPi =Ц
i=1
n-1
Е p+nPn|(1 -р)
i=1
Моменты искомого распределения можно получить многократным численным дифференцированием в нуле таблично заданной функции Б(в), в = 0, Н, 2Н, ... с последующей сменой знаков у нечетных производных. Однако предпочтительнее вычислять их непосредственно. Обозначим Е(х) набор моментов показательного распределения с параметром х, имеющий компоненты
Е1 (х) = I!/х1, I = 1, 2,...
Символом * будем обозначать свертку распределений в моментах. Для наборов моментов распределений между выходящими заявками сохраним индексацию прежних обозначений, заменив лишь 5(в) на й. В модели М/М/п при числе заявок в системе у > п
йу = Е(пц), у = п, п +1, ... .
Перепишем формулу для 8У (в) при у < п в виде
8j (s) =
І> + А
j' jjj, + A + s
i> № + X
ІМ
■1 + -
X
j\x + X jix + X
8j+l(s)
X № + Х г , .
-----------------т---------—------5 ■+1(в).
уц + А і^ + Х + в № + Х іц + Х + в ‘
Заменяя преобразование Лапласа показательных плотностей их моментами и произведение ПЛС - сверткой, убеждаемся, что
d; = j м+х
j E(& + Х) + -^E(jix + X)*dj,
в частности:
iV + A j = n -1, n - 2, ..., 0,
d0 = E(X)* d1.
Заметим, что свертка в моментах распределений А и В может быть выполнена на основе символического разложения
сп = (а + Ь)п
(9)
путем перевода показателей степени в индексы соответствующих моментов. Для выполнения сверток полезно составить автономную процедуру.
Результаты практически совпали с моментами распределения интервалов между входящими заявками, а коэффициенты немарковости с точностью до пяти знаков оказались нулевыми (выходящий поток - простейший). Это согласуется с классическим результатом Берке, позволяет считать основную идею расчета правильной и распространить ее на более сложные немарковские системы, для которых теория расчета выходящего потока до сих пор отсутствовала (в частности, на интересующую нас «дважды коксову»).
Коксова система. Введем обозначения для у-го яруса:
М,л - суммарная интенсивность обслуживания в £-м_микросостоянии каждой ветви;
М 1,1 - интенсивность переходов вверх по завершению обслуживания;
М ц - интенсивность переходов вправо по завершению фазы обслуживания (сумма двух последних интенсивностей равна Му,;);
С-1/, С(2) - наборы моментов распределения времени до ближайшего обслуживания при соответствующем исходном состоянии (верхние индексы указывают левую и правую ветвь диаграммы соответственно);
ау - количество микросостояний в одной ветви диаграммы.
В этих обозначениях для у > п и правой ветви диаграммы имеем
-М п,і
Мп,і + ^2
Е(Мп,і + ^2) +
Мп,і
-сі2? ■
-сі1)
МпЛ +^2 п,і Мп. +^2 п,і
і °п, ^п 1 • • •, 1 •
Для тех же значений у и левой ветви диаграммы надо дополнительно учитывать возможные переходы в правую ветвь по завершению фазы прибытия заявки:
Щ = .Мп,і, Е(Мп,і +Х1) + Е(Мп,і +Х1)*
Мп,і +Л1
Мпл _с(1) , и^1 с(2) | и^1 с(1)
Мпі +^1 п,і Мп. +^1 п,і Мп. +^1
і — Оп, Оп ~ 1 • • •, 1 •
После выполнения сверток в моментах согласно формуле (9) замечаем, что для каждого I искомые {Си,/] входят в обе части приведенных уравнений, причем «с перекрестом». Поэтому для их определения приходится последовательно решать ап систем из шести линейных алгебраических уравнений (по три момента в обеих ветвях на каждое микросостояние).
Для «ненасыщенных» ярусов у = п-1, п-2, ..., 1 в правой ветви диаграммы прибытие заявки приводит в одноименные состояния левой ветви нижележащего яруса. Соответственно имеем рекуррентные формулы
с(2> = Му,і Е(Му і + х2) +
у,і Муі +^2 ( у,і 2)
+Е(Мц +Х2)*
Муі -с<2> +-
Хо
2 с(1)
М:і +х9"у,і+1т м;і +х. 1+1,і
у,і ' ' >2 1иі,і^' 12
і — О у, О у ~1, • •., 1 •
Для левой ветви зоны ненасыщенных состояний завершение фазы прибытия заявки может привести как в правую ветвь диаграммы, так и в нижележащий ярус:
У,і Міл +*1
Е(Мц +Х1) + Е(Мц +Х1
Му,і _с« + иЯ1 с<2} + иЯ1 с£>
Мц +Х1 у,і+1 ■ Му і +ЛГу,і ■ Му ^ +І1
І+1,і
і ау, ау 1, •••, 1
Здесь первое слагаемое соответствует непосредственному завершению обслуживания, а последующие - переходам по завершению первой фазы обслуживания и фазы прибытия заявки (окончательной или с переходом во вторую фазу ожидания прибытия) с рекуррентным учетом моментов распределения времени до ближайшего обслуживания во вновь возникшем состоянии.
Наконец, в случае полного освобождения системы до очередного завершения обслуживания сначала придется дождаться прибытия хотя бы одной заявки. Здесь
<1 = Е(^)* СЦ); С$ = Е(Х1)* (иС,*2! + йС^1!) с111) .
В «ненасыщенной» зоне искомые наборы моментов определяются последовательно через ранее вычисленные: для каждого I сначала вычисляются «правые», а затем «левые» моменты.
Распределение времени ожидания
Интересующее нас распределение можно задать его первым моментом уЬ1 и набором коэффициентов немарковости. Тогда высшие моменты можно вычислить по формуле
£ = £>1 (^ +1!), I = 2, 3, ... . (10)
Известно, что условное распределение времени ожидания заявки в системе О1/М/п подчинено показательному закону. Это обстоятельство позволяет считать, что вид распределения времени ожидания не зависит от распределения интервалов между смежными заявками и числа каналов и определяется только распределением времени обслуживания. В работе [4] предложен конструктивный метод решения этой задачи, основанный на сохранении коэффициентов немарковости системы М/О/1.
Проблемы программной реализации
Зависимость структуры матриц переходов от типа и порядка аппроксимирующих фазовых распределений ставит ценность программной реализации расчетной схемы в прямую зависимость от возможности автоматического построения матриц переходов. Эту проблему в общем случае можно решить следующим образом: ____
- для каждого у-го яруса системы, у = 0, п, автоматически сгенерировать последовательность ключей микросостояний;
- формировать матрицы переходов, сопоставляя «ключи-источники» у-го яруса и совместимые с ними по выбранному типу переходов ключи-«пре-емники» для матрицы Ву на (у-1)-м ярусе, для ¥ у -на у-м, для Ау - на (у + 1)-м.
В обсуждаемом случае С2-аппроксимации структура матриц и правила их формирования относительно очевидны. Однако программная реализация метода оказалась значительно сложнее, чем это представляется при ознакомлении с его идеей. Дело в том, что метод предполагает работу с переменным числом матриц перехода изменяемой размерности, определяемой в процессе счета. Соответственно при программировании на Фортране пришлось упаковывать матрицы каждого вида (и их половинки для треугольных матриц) в линейные массивы, выделять под эти объекты динамическую ([а11оеа1аЪ1е]) память, фиксировать отдельно их размеры и координаты в линейных структурах, явно переадресовывать операнды матричных операций.
На эти проблемы наложились порождаемые специфической структурой упомянутых матриц: А - блочная с ненулевыми блоками в виде диагональных матриц, В - блочно-диагональная с одинаковыми блоками, ¥ - верхняя треугольная, Б -диагональная.
Учет этих особенностей при программировании требуемых матричных операций оказался нетривиальным и весьма поучительным.
Тестирование метода
Для проверки корректности метода и его программной реализации оценивалась средняя длина очереди (показатель компактный, наглядный, чувствительный к входным данным и инвариантный к масштабу времени) при коэффициенте загрузки системы р = 0,7 для различных видов исходных распределений и числа каналов. Для моделей (в обозначениях Кендалла) М/М/1, М/Е3/1, М/Б/1, М/Г15/1, Е4/М/1, Б/М/1, Г15/М/1, М/М/2, М/Б/2, М/Ез,/2, Е4/Б/2 были получены соответственно значения 1,633, 1,089, 0,817, 1,361, 0,861, 0,599, 1,292, 1,345, 0,692, 0,911, 0,113. Эти результаты практически совпали с полученными ранее посредством других численных методов, в том числе не использующих фазовые аппроксимации (например, метода Кроммелина для систем М/Б/п, Быкадорова для Ек/Б/п, Такача для О1/М/п).
Правильность расчета выходящего потока подтверждается:
- равенством средних интервалов между заявками входящего и выходящего потоков;
- получением практически нулевых коэффициентов немарковости выходящего потока для моделей вида М/М/п, обработанных по общей схеме;
- приближением этих коэффициентов к таковым для распределения времени обслуживания -в случае одноканальных систем при устремлении коэффициента загрузки к единице;
- совпадением результатов обсчета вышеперечисленных моделей с полученными посредством менее универсальной аппроксимации Н2/Н2/п. Частные случаи и обобщения
Самостоятельную ценность имеет наиболее типичный частный случай рассмотренной модели - с простейшим входящим потоком. В этом случае диаграмма переходов сводится к одной ветви и программная реализация сильно упрощается.
Как обобщение алгоритма можно рассматривать его вариант для системы с ограниченной очередью. В этом случае в микросостояния нижнего яруса диаграммы нет переходов снизу. Соответственно для этого яруса корректируется матрица Б-¥ и отпадает надобность в расчете одного из вспомогательных векторов итерационного метода (Р*). Для определения вероятности свободного состояния на заключительном этапе итерационного метода здесь условие баланса заявок заменяется условием нормировки вероятностей. Дальнейшим обобщением является расчет замкнутой системы с интенсивностью входящего (простейшего) потока, зависящей от числа заявок в системе - в типичном случае линейно убывающей до нуля. Здесь дополнительно модифицируются числовые значения матриц.
Заключение
Конкретные результаты данной статьи сводятся к следующим.
1. Предложен алгоритм расчета параметров распределения Кокса С2, выравнивающий три заданных момента.
2. Проведено исследование алгоритма; выявлена область входных данных, порождающая «патологические» вероятности переходов (комплексные и отрицательные); показано, что при обратном переходе от вероятностей микросостояний к укрупненным вероятностям патология исчезает и конечные результаты согласуются с полученными существенно иными методами. Указаны преимущества С2-аппроксимации в сравнении с другими фазовыми распределениями.
3. Впервые предложены диаграммы переходов для «дважды коксовой» системы, на основе которой получены матрицы интенсивностей переходов.
4. Разработан способ расчета распределения интервалов между заявками выходящего потока,
Литература
1. Рыжиков Ю. И. Теория очередей и управление запасами: Учебник. СПб.: Питер, 2001. 384 с.
2. Рыжиков Ю. И. Алгоритм расчета многоканальной системы с эрланговским обслуживанием //Автоматика и телемеханика. 1980. № 5. С. 30-37.
3. Рыжиков Ю. И., Хомоненко А. Д. Итеративный метод расчета многоканальных систем с произвольным распределением времени обслуживания // Проблемы управления и теории информации. 1980. № 3. С.32-38.
4. Рыжиков Ю. И. Три метода расчета временных характеристик разомкнутых систем массового обслу-
в частности - непосредственно моментов этого распределения.
5. Даны рекомендации по программной реализации основных этапов алгоритма.
Эти результаты вместе с методами решения уравнений баланса и расчета временных характеристик позволяют, наконец, сравнительно просто и достаточно точно проводить полный расчет многоканальных систем с произвольными распределениями интервалов между заявками и времени обслуживания.
В сравнении с имитационным моделированием метод существенно выигрывает по объему и точности вычислений. Он избавляет от изучения и использования других фазовых аппроксимаций и (если не говорить о подготовке будущих теоретиков) -вообще всех прочих методов расчета немарковских систем обслуживания. Метод может использоваться как составная часть процедур расчета сетей обслуживания методом потокоэквивалентной декомпозиции [1].
живания || Автоматика и телемеханика. 1993. № 2. С. 127-133.
б. Evans R. D. Geometric Distribution in Some Two Dimensional Queuing Systems || Operat. Res. 19б7. Vol. 1б. N б. P. 830-84б.
6. Khomonenko A. D., Bubnov V. P. A Use of Coxian Distribution for Iterative Solution of M|G|n|R < <»\le \infty queuing systems || Problems of Control and Information Theory. 198б. Vol. 14. N 2. P. 143-1б3.
7. 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. 197б. Vol. 19. N 2. P. 147-1б7.
