Итерационный метод расчета системы с коксовским обслуживанием M/C2/n и его программная реализация Текст научной статьи по специальности «Математика»

rimental verification of the proposed identification method are presented, showing the possibilities of its use in the technological process of information objects automated control, as well as its effectiveness in relation to existing methods.

Key words: identification, binary data arrays, certification tests, big data.

Lebedenko Evgeniy Viktorovich, candidate of technical sciences, docent, lebeden-ko_eugene@mail. ru, Russia, Orel, Academy of the Federal Guard Service of Russian Federation,

Ryabokon Vladimir Vladimirovich, employee, mimicria@,mail. ru, Russia, Orel, Academy of the Federal Guard Service of Russian Federation,

Lapko Aleksandr Nikolaevich candidate of technical sciences, employee, lan46@,mail.ru, Russia, Orel, Academy of Federal Guard Service of the Russian Federation,

Kutsakin Maksim Alekseevich, employee, max_kooks@,mail. ru, Russia, Orel, Academy of the Federal Guard Service of Russian Federation

УДК 519.872

ИТЕРАЦИОННЫЙ МЕТОД РАСЧЕТА СИСТЕМЫ С КОКСОВСКИМ ОБСЛУЖИВАНИЕМ M/C2/n И ЕГО ПРОГРАММНАЯ РЕАЛИЗАЦИЯ

А.В. Уланов, В. А. Лохвицкий, Д.Ю. Старобинец, А.Ф. Шинкаренко

Рассматривается система массового обслуживания с коксовским обслуживанием MIC2ln. Представлена технология построения графа марковизированной системы, матриц интенсивностей инфинитезимальных переходов и алгоритм расчета стационарного распределения числа заявок. Рассмотрены особенности хранения матриц и реализации вычислений на языке Фортран. Показана возможность применения распределения Кокса второго порядка с комплексными и парадоксальными параметрами. Результаты расчета верифицированы с помощью имитационного моделирования.

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

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

Эффективным методом марковизации многоканальных систем массового обслуживания (СМО) является применение распределений фазового типа (обозначаются РИ). Упомянутые распределения представляют собой параллельно-последовательные комбинации фаз прохождения заявок с экспоненциально распределенной задержкой в каждой из них. Идея метода фиктивных фаз, когда исходные неэкспоненциальные распределения заменяются фазовыми, была выдвинута еще А.К. Эрлангом. Порядок аппроксимации при этом может определяться количеством сохраненных начальных моментов исходного распределения.

Наиболее распространённые распределения фазового типа - Кокса, Эрланга и гиперэкспоненциальное. Применение того или иного распределения зависит от вероятностной картины процессов поступления и обслуживания заявок. Так, для аппроксимации распределений с коэффициентом вариации и > 1, как правило, используют гиперэкспоненциальную (И) аппроксимацию, а в остальных случаях - эрланговскую (Е). При этом предполагается, что параметры распределений принимают только вещественные значения. Тем не менее, в работах [1, 2] показана применимость гипр-экспоненциального распределения второго порядка с комплексными параметрами для аппроксимации времени обслуживания и других входных характеристик СМО с произвольным (в том числе меньшим единицы) коэффициентом вариации.

В представленной статье показана возможность комплексных параметров для распределения Кокса второго порядка (С2), обсуждается технология построения диаграммы переходов марковизированной СМО М1С21п, особенности формирования матриц интенсивностей инфинитезимальных переходов и программной реализации итерационного алгоритма расчета стационарного распределения числа заявок.

1. Цепь Маркова для системы М1С21п

При коксовском распределении второго порядка заявка, попадая в канал обслуживания, обязательно задерживается в первой экспоненциальной фазе с интенсивностью ц1, а затем с вероятностью у попадает во вторую экспоненциальную фазу с интенсивностью обслуживания или с вероятностью у = 1 - у покидает систему. На рис. 1 представлена схема этого распределения.

Рис. 1. Распределение Кокса второго порядка

Поскольку С2-распределение имеет три параметра, оно позволяет выровнять три начальных момента исходного, чего для инженерных расчётов этого вполне достаточно.

Параметры С2-распределения легко получить аналитически. В соответствии со схемой (рис. 1) можно вывести уравнения для вычисления начальных моментов:

1 У /■

—+ —= /ь

т т2

1 /1 1ч/-

"У + У (-+ = f2,

т2 т1т2 т2

1/1 1 1ч/-

-Т + +-2 + _3) = /з.

_ т 2 №2 т2

Здесь gi = /¡/И, I = 1,3 - начальные моменты распределения Кокса второго порядка.

Параметры С2-распределения находятся аналитически:

„ = (/1/2 - /3 )

т 2 =-2-,

2(/22 -/1/3)

т =//, (1)

т 2 /2- /1

= т 2(/т -1) т ,

где

V = (/1 /2 - /з)2 - 4(/22 - /1 /3)(/12 - /2).

Очевидно, что в случае V < 0 параметры С2-распределения принимают комплексные значения.

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

На рис. 2 приведён пример диаграммы переходов для трехканаль-ной СМО М/С2/3.

Столбец справа на диаграмме - число заявок в системе. Двухразрядные числа соответствуют расстановкам заявок по фазам С2-обслуживания и называются ключами микросостояний. В соответствии со схемой распределения (рис. 1) первый разряд - число заявок, находящихся в первой фазе, второй разряд - число заявок во второй фазе. Сумма цифр каждого ключа микросостояний всегда равна числу заявок в каналах обслуживания.

Рассмотрим всевозможные переходы из микросостояния с ключом «11» на втором ярусе. Значение ключа говорит о том, что в нашей трехка-нальной системе две заявки, в одном канале заявка первой фазе С2-

279

обслуживания, в другом - во второй фазе. Вновь прибывающая заявка попадает в свободный канал обслуживания в первую фазу, поэтому с интенсивностью X система перейдет в микросостояние с ключом «21». С интенсивностью заявка из первой фазы с вероятностью у покинет канал, число заявок в системе уменьшится на одну, а система перейдет на вышележащий ярус в микросостояние с ключом «01». С вероятностью у заявка останется в канале, но перейдет во вторую фазу, в результате чего значение ключа станет равным «02». Наконец, с интенсивностью заявка из второй фазы покинет канал и система перейдет в микросостояние с ключом «10».

При полной занятости диаграмма стабилизируется на всех нижележащих ярусах.

В зависимости от значений выравниваемых моментов параметры С2-распределения могут принимать комплексные и «парадоксальные» значения. В статье Ю.И. Рыжикова [3] приведено исследование С2-аппроксимации исходного гамма-распределения с коэффициентом вариации и, было выявлено, что:

- случай 0 < и < 1/ 42 дает комплексные параметры;

- при 1 / л/22 < и < 1 параметры вещественные;

- при и > 1 имеем параметры С2-аппроксимации вещественны, но парадоксальны: «вероятность» у прохождения двух фаз будет меньше нуля.

Поскольку параметр у С2-распределения обслуживания интерпретируется как вероятность ухода заявки из канала после прохождения первой фазы обслуживания, большинство специалистов по теории массового обслуживания рассматривают лишь тот случай, когда данный параметр определен на вещественном интервале [0; 1], что соответствует аппроксимации распределений с коэффициентом вариации 1 / л/2 < и < 1.

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

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

2. Решение уравнений баланса переходов

Обозначим Б' множество всех микросостояний системы на /-м ярусе, когда на обслуживании находится ровно / заявок, а через о/ - количество элементов в Б/. Далее на основе диаграммы переходов определим следующие матрицы переходов инфинитезимальных переходов:

А'[°/'х°/'+Л - в 8/+1 (прибытие заявки),

В[о/хо/.;] - в Б,.! (полное завершение обслуживания заявки),

С/'[а/ха/'] - в Б,- (конец промежуточной фазы обслуживания),

Ц'[о/ХО/] - ухода из состояний яруса/ (диагональная матрица).

При/ > п матрицы перестают меняться ( п - число каналов).

Введём векторы-строки вероятностей у/ = (у/1,у/1,...,у/ о/ } нахождения СМО в стоянии (/', /),/ = 0, 1,... Теперь можно записать вектор-но-матричные уравнения баланса переходов между микросостояниями диаграммы [1]

7 0 ^о = 71В l,

— (2) = 7/-1А-1 + 7/С/ + 7/+1В/+1, ' =1, Я

где Я - число обсчитываемых ярусов диаграммы.

Для системы М/С2/3 матрицы инфинитезимальных интенсивностей выглядят следующим образом:

=А';

"1 "1 о о о

о о о"

А II "1 о о" 1 о 1 о о

, А2 = о о о А <-о II

о 1 о_ 1 о о 1 о

о о о о о о 1

В1

¡¡1 у

.т 2.

В

С1

о ¡¡1 у о о

, С 2

2т1у о ' ¡¡2 т2У о 2^2

о 2^1у о

В

о о

¡¡1 у о

'3[цу о о

¡¡2 2^1у о

о 2т2 т1у

о о зт2

о з^у о

, с3

В

о

о о 2т1у о о о о ¡цу о о о о

С/;

Б = [1], £>1 =

1 + т1

о

о

1 + т 2

в2 =

1 + 2^1

0

о 1 + ^1 +т 2

0

о

б =

1+зт1 о о о

о

о

о

1+2т1 +т 2 о о

о 1+т1 + 2т 2 о о о 1+зт2

=

о о

1 + 2т 2

у > з.

Одним из достоинств распределения Кокса является простота выявления закономерностей для автоматического построения матриц интен-сивностей.

Система уравнений (2), дополненная условием нормировки, для многоканальных СМО характеризуется чрезвычайно большой размерностью, и стандартные методы решения систем линейных алгебраических уравнений применительно к ней оказываются малоэффективными. Для ее решения применяют метод матрично-геометрической прогрессии, впервые предложенный Ивэнсом [4] и развитый школой Ньютса [5], либо итерационный метод, предложенный Такахаси и Таками [6]. Остановимся подробнее на последнем.

Положим, что г у = уу/ру , где ру - суммарная вероятность наличия в

системе ровно у заявок, и обозначим

Ху = Ру+11 Ру , = Ру-11 Ру .

Тогда систему (1) можно переписать относительно векторов условных вероятностей {гу}, нормированных к единице в пределах яруса:

г оБо = г оСо + x0t1B1, (3)

= 2угу-1 Ау-1 + г]С] + хугу+1ву+1,

у = 1,2,

т

С помощью векторов-столбцов 1у = {1,1,.. .,1} размера су для всех у дополним систему (3) уравнениями нормировки

Ч1у =1 (4)

и баланса суммарных интенсивностей переходов между смежными ярусами

гуАу 1у+1 = хугу+1 ву+11 у.

(5)

Алгоритм расчёта набора векторов {гу} и чисел {ху} и {¿у}, удовлетворяющих соотношениям (3)-(5), в случае разомкнутой системы с неограниченной очередью опирается на существование предельного при у ^ да вектора условных вероятностей гда, что является следствием стабилизации матриц интенсивностей при у > п. Алгоритм основан на последовательном

приближении к искомым характеристикам для ограниченного множества

индексов у = 1, Я и по существу является блочным вариантом известного

метода Гаусса-Зейделя.

Перепишем уравнения системы (3) для у > 1 в виде

гР(Оу - Су) = ¿("У1» А,-! + х<"> г" в+ъ у = 1,2,..., где верхний индекс показывает номер итерации. Тогда можем записать

г ( ) = 7

" Р) + х " Р} , (6)

где

Р у = г т-1) А-1(О - Су)Л (7)

Р} =-Су)-1.

При у = Я считается, что

= гЯ"-1) Вя (ОЯ - Ся)-1. (8)

Для расчета {х(т)} и {г(т)} перепишем (5) с учётом (6):

(7(")Р' + х(")Р" ) А 1 = х(")г (") В 1

(2] Р}+ х] Р}) А7 17+1 = х] ]+1В у+1 у. Отсюда следует пропорциональность

с коэффициентом

z (т) = сх(т) (9)

с ' j+1 BJ^J -Р"A 1J+ (10)

p J aj1J+1 '

В этой и последующих формулах произведения матриц переходов на вектора 1j равны суммам строк и вычисляются до начала итераций. Подстановка (9) в (6) и умножение обеих частей на 1j дают

1 = t fl J = х f(Cp J +Р" )1 ,.

Итак,

х m=i/[(cp j+p" )i j ]. (ii)

Условие прекращения итераций

(12)

х (m) - х (m-1) xj xj

тах

7

В качестве начальных приближений целесообразно задать Я = г¥. тогда существуют предельные значения отношений смежных вероятностей х и г, причём из (5) следует

х = гА1п/гВ1п 283

и стационарные вероятности с большими индексами образуют геометрическую прогрессию. Расчёты свидетельствуют о хорошей аппроксимации x формулой

х = р2/(иА ),

где р - коэффициент загрузки системы; иА, иВ - коэффициенты вариации интервалов между заявками входящего потока и времени обслуживания соответственно.

Максимальное число вероятностей Я находится по формуле

Я » п + (1§2 - s), где п - число каналов. При этом учитываются состояния, вероятность которых не превышает Ю-8.

После прекращения итераций можно переходить к нахождению стационарных вероятностей. Из определения (х/} следуют равенства

Р/+1 = р/х/, (13)

Из закона сохранения заявок

п-1

Е (п - /)Р/ = п - Ь/а,

/=о

получаем формулу для вероятности свободного состояния системы

п - Ь а

Р о = п-1 /-1 . (14)

п + Е(п - ') П х /=1 /=о

Последующие вероятности для / = 1, Я определяются рекуррентно с помощью (13). При необходимости та же формула может быть применена для больших значений / с использованием х/ = хда.

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

3. Организация структур данных

Для хранения матриц интенсивностей, векторов вероятностей микросостояний и других массивов переменной длины выделялась динамическая (ЛЬЬОСЛТЛВЬЕ) память. Матрицы интенсивностей определялись как трёхмерные массивы, вектора вероятностей микросостояний - как двумерные.

Матрицы А' интенсивностей прибытия заявок являются диагональными и их умножение на другие матрицы сводится к умножению каждого элемента на интенсивность входящего потока. Поэтому в выделении под них памяти необходимости нет.

Под прочие матрицы переходов в процессе вычислений выделялись массивы B(n,n+1,n+1), C(0:n,n+1,n+1), D(0:n,n+1,n+1) с комплексным типом данных complex*16.

При обращении матриц (Df - Cj) необходимо было учесть, что они являются верхними треугольными.

Матрицы Bj + 1(Dj - Cj), j = 1, n рассчитывались до начала итераций и заносились в массив BD(0:n,n+1,n+1).

Матрицы Bj _ 1(Dj - Cj), j = 1, n+1 рассчитывались до начала итераций и заносились в массив AD(n+1,n+1,n+1) 4. Программная реализация

Расчет стационарного распределения числа заявок в системе производился в процедуре cox2(lam,mu,y,n,jmax,eps,p). Исходными данными для процедуры расчета стационарного распределения числа заявок в системе являлись:

_ интенсивность входящего потока lam, _ параметры распределения Кокса mu(2) и у,

_ требования к точности вычисления отношений смежных вероятностей eps,

_ количество вычисляемых вероятностей jmax, _ число каналов n.

После прекращения итераций вычисляется распределение числа заявок p, передаваемое основной программе.

В качестве среды программирования выбрана система Simply Fortran, которая работает под управлением Windows 7.

Основными этапами расчета являются генерация матриц интенсив-ностей, расчет начальных приближений векторов {j и {х}, их уточнение и расчет стационарного распределения числа заявок {pj}.

На начальном этапе происходит генерация матриц интенсивностей. В качестве примера приведем фрагмент кода процедуры cox2, генерирующий матрицы В{:

do i=1,n k=i do j=1,i b(ijj)=k*mu(1)*(1-y) b(i,j+1,j)=j*mu(2) k=k-1

end do end do.

Внешний цикл задает номер яруса, для которого генерируется матрица, внутренний _ выполняет непосредственно генерацию.

После генерации следует расчет начальных приближений векторов условных вероятностей микросостояний, хранящихся в массиве t(0:jmax,n+1). Возможно задать начальные приближения равновероятными.

Далее следует цикл с условием do while (delta>eps), контролирующий выполнение условия (12), в который вложен цикл с параметром по номеру ярусаJ, в котором осуществляется уточнение элементов xj и j Основными этапами при этом являются:

I к

- расчет векторов Р j и Р j по формуле (7),

- вычисление коэффициента с по формуле (10),

- расчет xJ по формуле (11),

- расчет z по формуле (9),

- уточнение tf по формуле (6),

- пересчет невязки в модуле выражения (12).

После прекращения итераций производится расчет стационарного распределения {pj} по формулам (13)-(14). Текст программы: xl=1d0 ' г=0

do j=0,jmax P(j)=xl xl=xl*x(j) end do m=0 do j=0,n-1

m=m+(n-j)*p(j) end do m=(n-lam*f(1))/m r=0

do j=0,jmax p(j)=m*x1(j) r=r+x1(j) end do. Здесь f(1) - среднее время обслуживания. 5. Численные результаты

Выполним расчет системы M/G/n на имитационной модели и численным методом через С2-аппроксимацию для следующего набора исходных данных: число каналов обслуживания n = 8, средняя интенсивность входящего потока X = 6.4, среднее время обслуживания заявок b1 = 1. Время обслуживания представлено следующими законами:

- вырожденным D (коэффициент вариации и = 0);

- равномерным U на интервале [0; 2] (и ~ 0,577);

- экспоненциальным M (и = 1);

- гамма Г с параметром формы 0,5 (и ~ 1,41);

- Вейбулла W c параметром формы 0,46 (и = 2,5).

В программе, реализующей вычисление распределения числа заявок итерационным методом, параметры С2-распределения подбирались по трем моментам упомянутых распределений длительности обслуживания -формула (1). Параметры С2-распределения, рассчитанные по трем начальным моментам исходного распределения B(t), приведены в таблице.

Параметры С2-аппроксимации

Б(0 У Ц1 ^2 Б(0 У Ц1 ^2

В 1,667+ /•0,471 2,0-/•1,414 2,0+/-1,414 Г0,5 0,414 3,414 0,586

и 1,0 1,5-/•0,866 1,5+/-0,866 Ж 0,238 3,673 0,327

Результаты расчета распределения числа заявок в системе для В-, и-, М-, Г- и Ж-распределений приведены на рис. 3.

Пунктирными линиями показаны графики, полученные с помощью имитационной модели, сплошными линиями - с помощью итерационного метода через С2-аппроксимацию.

Рис. 3. Распределение числа заявок в системе М/С2/8

Из графиков видно согласие результатов даже в области комплексных и параметров С2-распределения. Расстояние Колмогорова между результатами, полученными методом имитационного моделирования и с помощью С2-аппроксимации В-, и-, М-, Г- и Ж-распределений времени обслуживания составило {0,008; 0,05; 0,003; 0,04; 0,004} соответственно, что говорит о приемлемой точности аппроксимации. Заметим, что для В- и и-распределений параметры С2-распределения принимали комплексные значения, что, тем не менее, не повлияло на корректность конечного результата.

Заключение

Распределение Кокса второго порядка может с успехом применяться для итерационного расчета СМО наряду с гиперэкспоненциальным, позволяет выровнять три начальных момента. При этом возникающие в некоторых случаях комплексные параметры не влияют на корректность конечного результата - распределение числа заявок.

Представленные в статье особенности программной реализации расчета системы M/C2/n могут применяться в учебных целях при изучении структур и алгоритмов обработки данных, поскольку диаграмма переходов марковизированной СМО является взвешенным ориентированным графом, а технология построения матриц интенсивностей, их хранение и обработка раскрывает особенности обработки графов.

Список литературы

1. Рыжиков Ю.И., Уланов А.В. Применение гиперэкспоненциальной аппроксимации в задачах расчета немарковских систем массового обслуживания // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2016. № 3 (36). С. 60-65.

2. Лохвицкий В. А., Уланов А.В. Численный анализ системы обслуживания с гиперэкспоненциальным «охлаждением»// Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2016. № 4(37). C. 36-43.

3. Рыжиков Ю.И. Полный расчет системы обслуживания с распределениями Кокса // Информационно-управляющие системы. 2006. № 2. C. 38-46.

4. Evans R.D. Geometric Distribution in Some Two Dimensional Queuing Systems // Operation Research. 1967. V. 15. No. 5. P. 830-846.

5. Neuts M.F. Matrix-Geometric Solutions in Stochastic Models: an Algorithmic Approach. Baltimore: J. Hopkins Univ. press, 1981.

6. Takahashi Y., Takami Y. A Numerical Method for the Steady-State Probabilities of a GI/G/c Queuing System in General Class // Journal of the Operations Research Society of Japan. 1976. V. 19. No. 2. P. 147-155.

Уланов Александр Викторович, канд. техн. наук, старший научный сотрудник, ulanov246@rambler. ru, Россия, Санкт-Петербург, Военно-космическая академия имени А. Ф.Можайского,

Лохвицкий Владимир Александрович, канд. техн. наук, докторант, vovan296@mail. ru, Россия, Санкт-Петербург, Военно-космическая академия имени А. Ф.Можайского,

Старобинец Дмитрий Юрьевич, преподаватель, norfl2arambler.ru, Россия, Санкт-Петербург, Военно-космическая академия имени А. Ф.Можайского,

288

Шинкаренко Антон Федорович, канд. техн. наук, начальник научно-исследовательской лаборатории, tonio8 7amail. ru, Россия, Санкт-Петербург, Военно-космическая академия имени А. Ф. Можайского

AN ITERATIVE METHOD OF M/C2/n QUEUING SYSTEM CALCULATING AND ONES PROGRAM REALIZATION

A.V. Ulanov, V.A. Lokhvitsky, D.Yu. Starobinets, A.F. Shinkarenko

A queuing system M/C2/n is shown in article. The Markow chain for this queuing system, technology of creating transmission matrixes and numerical iterative method for calculating of stationary probabilities of number of customers are presented. In this article is described program realization of numerical method for FORTRAN. The possibilities of Coxian distribution of second order with complex parameters are shown. The results are verified using simulation.

Key words: numerical method, Coxian distribution, program realization, non-markovian queuing systems.

Ulanov Alexander Viktorovich, candidate of technical sciences, senior research associate, [email protected], Russia, St.-Peterburg, Mozhaisky Military Space Academy,

Lokhvitsky Vladimir Aleksandrovich, candidate of technical science, doctoral candidate, vovan296amail. ru, Russia, St.-Peterburg, Mozhaisky Military Space Academy,

Starobinets Dmitriy Yurievich, teacher of computer software department, norfl 2 aram bler. ru, Russia, St.-Peterburg, Mozhaisky Military Space Academy,

Shinkarenko Anton Fedorovich, candidate of technical science, chief of research laboratory, tonio87a mail.ru, Russia, St.-Peterburg, Mozhaisky Military Space Academy