CC BY
83
УДК 519.248, 519.713.2, 519.178
КОМПЛЕКС ПРОГРАММ ИМИТАЦИОННОГО МОДЕЛИРОВАНИЯ РАБОТЫ СИСТЕМЫ МАССОВОГО ОБСЛУЖИВАНИЯ С НЕОДНОРОДНЫМИ ПРИБОРАМИ И РАЗДЕЛЬНЫМИ НАКОПИТЕЛЯМИ
А. П. Котенко, М. Б. Букаренко
Самарский государственный технический университет,
Россия, 443100, Самара, ул. Молодогвардейская, 244.
E-mails: akol959@mail. ru, maxim. bukarenko@gmail. com
Представлен комплекс программ имитационного моделирования работы и автоматического построения графа состояний системы массового обслуживания с неоднородными приборами и раздельными накопителями. В основе алгоритма лежит разработанный авторами метод аналитического описания таких систем конечными автоматами Мили.
Ключевые слова: система массового обслуживания, неоднородные приборы, имитационное моделирование, распределение интервалов времени между заявками, раздельные накопители, визуализация графа, матрица смежности графа.
Введение. В нотификации [1] рассмотрим систему массового обслуживания (СМО) с к > 0 неоднородными приборами и раздельными очередями сигнатуры Т = Т(т\ х /XI; т 2 х Ц2', ■ ■ ■] тпк х №к), где — пропускная способность, т,г - число мест в очереди г-того прибора, г € 1 ,к. Её граф уже не будет графом процесса гибели-размножения, но будет определяться протоколом диспетчеризации входных заявок между приборами.
В распространенных на практике случаях аналитические методы моделирования данных СМО не работают.
В первую очередь, это связано с тем, что они распространяются только на системы с простейшими потоками событий. Уравнения Колмогорова для марковских процессов составляются при условии существования предела
ду = цт Ш
13 д*->0 А* ’
где рц — вероятность перехода системы из г-того в ^’-тое состояние. Это влечёт за собой пуассоновское распределение интервалов времени между событиями потока.
Кроме того, моделирование СМО с использованием тем или иным образом уравнений Колмогорова подразумевает отсутствие последействия, когда для любых двух непересекающихся отрезков времени число событий, попавших на один из них, зависит только от длины отрезка, но не зависит от числа событий, попавших на другой отрезок. Функциональная связь между моментами появления новых заявок в СМО с самым простым регулятором нарушает данное условие.
Андрей Петрович Котенко (к.ф.-м.н., доц.), доцент, каф. прикладной математики и информатики. Максим Борисович Букаренко, аспирант, каф. прикладной математики и информатики.
При этом уравнения Колмогорова эффективно решаются только в случае однородных марковских цепей, когда гарантированно существует стационарное решение. Наконец, немаловажна громоздкость системы уравнений Колмогорова для СМО с различимыми каналами, порождённая большим числом состояний системы.
Представление такой СМО конечным автоматом (КА) Мили позволило создать комплекс программ для моделирования работы системы с произвольным распределением интервалов времени между появлением входящих/обработанных заявок, который строит матрицу смежности и граф состояний данной СМО как без диспетчеризации (т. е. распределение входящих заявок по доступным приборам произвольно и не зависит от их загруженности), так и в соответствии со следующим протоколом диспетчеризации.
Пусть очередная заявка обнаруживает систему в состоянии (х\, Х2, ■ ■ ■,
2/1, у2, ■ ■ ■, Ук), не являющемся состоянием отказа (1,1,... ,1; т\, т,2,..., т&). Если существует единственный прибор (с номером г), способный принять заявку (0 ^ Уг ^ |ХгГГц — 1|), то заявка направляется к нему. В противном случае выбираем г-тый прибор, способный обработать заявку с минимальным средним суммарным временем Т обслуживания напрвленных к нему заявок:
Т= ШШ 11~1{уг + 1+Хг), (1)
г:05С;у;5С|11:;т,;-1|
где случайная величина Хг £ [0; 1] характеризует незавершённость обработки заявки г-тым прибором в момент поступления новой заявки.
Программный комплекс содержит следующие программы 1-3.
1. Программа имитационного моделирования потоков входящих и обработанных заявок. В основе работы программы лежит алгоритм [3,4].
При имитации реальной СМО искусственные реализации входящего и выходящего потоков необходимого объема традиционно предполагают лишь пуассоновское распределение интервалов времени между входящими и обработанными заявками. С помощью КА можно сохранить не только произвольные законы распределения, но и структуру входного/выходного сигнала, ПОД которым понимается эмпирический ЧИСЛОВОЙ ряд [Т\, ?2, . . . , Тг, . . . , Тд.] из интервалов времени Тг между г-той и (г — 1)-й по счёту буквой входного алфавита (1—для входящего потока, или 0 — для потока обработанных заявок) [2]. Применив к этому ряду методы цифровой обработки сигналов, получим оценку необходимого числа искусственных реализаций.
На первом этапе выделим характерные циклы моделируемых потоков входящих/обработанных заявок: например, для абонентской активности оператора сотовой связи периоды циклов будут составлять день, неделю, год и т.п. Для этого методом максимальной энтропии Берга оценим спектральную плотность мощности (СПМ) каждого эмпирического сигнала.
Выбор порядка р авторегрессионной модели производится по критерию длин минимального описания [3]:
МЬО[р] = ЛПпрр +р1пЖ,
где сц — параметры модели Берга, Ж —длина сигнала, рг = рг~ 1(1 — |ецМ|2) монотонно уменьшаются с ростом порядка модели.
Выбор частоты отсечки производится по впадинам графика СПМ так, чтобы в полосу пропускания попали все доминирующие циклы. Далее сигнал очищается медианной фильтрацией от импульсных шумов, когда перепады значений сигнала велики по сравнению с дисперсией аддитивного белого шума. Выходной временной ряд ук скользящего медианного фильтра с окном ширины 2п + 1 для /с-того отсчета формируется из входного эмпирического ряда (... ,хк-1,хк,хк+1, ■ ■ •) формулой
У к •— тес! (Хк—п, Хк—п-\-\, . . . , Хк—\, Хк, Хк-\-1, • • • , Хк-\-п— 1) %к-\-п) )
где тес! (х\,..., Х2п+1) := хп+\ при ранжировании значений попавшего в окно вариационного ряда:
Ж(1) = тт (х!,х2, • • •, х2п+1) ^ ж(2) ^ ^ Х(2п+1) = т&х (х1,х2, ■■■, х2п+г).
Если эмпирический сигнал является суммой = у^-\- ^ детерминированной и стохастической величин, то сохранение его структуры при имитационном моделировании означает представление каждой искусственной реализации суммой Т( = у г + где ^ — случайная величина с тем же спектром дисперсий, что и При этом в качестве тренда у г может выступать как отклик низкочастотного линейного фильтра (например, в соответствии с описанной выше процедурой медианной фильтрации), так и заданная аналитически функция (найденная, к примеру, методами асинхронного гармонического анализа).
Моделирование стохастической составляющей £(£) производится по спектральному разложению
к=О
при некоррелированности центрированных случайных величин ик, Ук, дисперсии которых для одинаковых индексов к составляют величины
1 [т 2 [т кттт
В° = т]0 к^г1т’ Вк = т]0 к*(т')с°8~т~с1т’
где /^ — автокорреляционная функция £(£) на (—Т;Т):
СЮ 7
к(т) = ^2ПкСоё^~-к=О
Данный метод позволяет моделировать работу СМО при наличии неслучайных составляющих в интенсивностях входного/выходного потоков, что нарушает условие отсутствия последействия. Программа создает требуемое число искусственных реализаций входящего потока заявок и потоков обработанных заявок для каждого отдельного прибора СМО.
2. Программа имитационного моделирования работы СМО. Эта часть программного комплекса реализована на языке Ма'ЫаЬ следующим образом.
Период моделирования разбивается на оптимальное с точки зрения соблюдения условия ординарности и удобства анализа результатов число тактов. В каждый такт вырабатывается п + 1 сигнал, где п — число приборов системы. Состояние системы определим вектором
*51ХП = {х\ + У1,Х2 + У2, ■ ■ ■ ,Хг + Уг, . . . ,ХП + уп} ^ Уг ^ ТП^, Xi € [0, 1],
в соответствии с принятой нотификацией. Каждый сигнал прихода/обработки заявки является двоичным числом «у. Описание сигналов представлено в таблице.
Описание сигналов имитационной модели
Порядковый номер сигнала в такте Значение Sij Описание
3 = 1 1 Появление входящей заявки. Если = т\ + 1, то2 + 1, ..т* + 1, ..тп + 1, то 54+1 = St, заявка получает отказ. В противном случае она направляется на свободный прибор в соответствии с принятой диспетчеризацией.
3 = 1 0 й+1 = я.
j = 2, п + 1 1 Выработка ;-тым прибором обработанной заявки: Га^+%-1,а^ + % > 1, 0+1 0+1 \0,ж^+%=0
j = 2, п + 1 0 й+1 = й
На основе разработанной нотификации и представления СМО конечным автоматом Мили на языке Ма'ЫаЬ была разработана система имитационного моделирования работы СМО с различимыми каналами для произвольных параметров потоков входящих и обработанных заявок.
Она позволяет учитывать следующие характеристики (рис. 1): распределение интервалов времени между появлениями входящих заявок; интенсивность потока входящих заявок (для пуассоновского распределения) и математическое ожидание интервала времени между входящими заявками (для прочих типов распределения); дисперсию потока входящих заявок; число приборов системы (произвольное конечное); распределение интервалов времени между обработанными заявками, покидающими систему; интенсивности потоков обработанных заявок (для пуассоновского распределения) и математические ожидания интервала времени между обработанными заявками (для прочих типов распределения) для каждого прибора системы; дисперсии потоков обработанных заявок для каждого прибора; емкости накопителей приборов (произвольные конечные).
Программа позволяет получить набор состояний исследуемой СМО на каждом такте в течение заданного пользователем времени работы имитационной модели. Все время моделирования Т необходимо так разбить на N равных тактов времени т, чтобы выполнялись следующие условия:Жт ^ Т; Р(п > 1 ,т) = 0. Второе условие означает, что вероятность появления более одной входящей заявки за один такт равна нулю, что обеспечивает ординарность потока входящих заявок. Аналогичным образом соблюдается ординарность и для потоков обработанных заявок.
Для выполнения первого условия сначала полагаем длину вектора генерируемых псевдослучайных чисел равной N\ = Т/М\ [г], затем увеличиваем значение N\ на единицу до тех пор, пока сумма сгенерированных псевдослучайных чисел не окажется больше времени работы модели Т: хг ^ Т. При
этом длина такта Tj = k min^ Zi, где к < 1 и задано пользователем.
Отметим, что минимально необходимый для ординарности потока промежуток времени между сигналами выработки заявки
т°гЧ = min{Tou*} г€[1,»г]
может быть меньше продолжительности такта тгп, определяемого исходя из характеристик потока входящих заявок. В этом случае полагаем т = min{тгп,тогЛ}.
Разработанная программа имитационного моделирования СМО с различимыми каналами позволяет найти следующие характеристики системы: динамика загруженности j-того прибора; доля отказов; динамика вероятности простоя.
Средняя загруженность j-того прибора рассчитывается как отношение математического ожидания числа заявок, обслуженных j-тым прибором за данный промежуток времени, к максимально возможному числу заявок, которое прибор в состоянии обработать за данный временной интервал:
(n+l)fci+i
L4, = („ + l гЧкш-к,)-' v
III 7 —Г~ 1
г=(»г+1)Ач
где kj — номер такта.
Доля отказов рассчитывается как отношение числа заявок, заставших систему в состоянии St = {mi + 1, m2 + 1,..., nii + 1,..., mn + 1} и пришедших
к ■ _L 7
за время j ^ = const, к суммарному числу заявок, пришедших в систему
за указанный интервал времени.
Пример динамики значений данных характеристик представлен на рис. 2.
3. Автоматическое построение графа состояний СМО. Разработанная программа построения графа состояний СМО решает следующие задачи: построение списка допустимых состояний; расчет и визуализация матрицы смежности вершин; расчет координат вершин для наиболее эффективной визуализации орграфа состояний СМО; визуализация графа.
Ее интерфейс реализован на VBA for MS Excel, а визуализация осуществляется специально разработанным приложением Graph Theory Toolbox для Matlab. Обмен данными между MS Excel и Matlab происходит посредством специального набора инструментов Excel Link Toolbox, позволяющего вызывать функции Matlab из модуля VBA for MS Excel.
Ввод параметров
- Mi
— Моделирование си тало в на приход заявки -Интенсивность, 1/M[t]
Сред не квадратическое отклонение
Время моделирования
Коэффициент длины такта, к
0.5
— Моделирование си тало в на выработку заявки Количество приборов Интенсивности, 1/M[t]
Среднеквадратические отклонения Емкости накопителей
[0.15 0.05 О С
[2 4 6]
Запуск
Рис. 1. Окно ввода параметров
1 Т = ГЗСИаи
Рис. 2. Результаты моделирования СМО сигнатуры Т = (2 х 0,15;4 х 0,05; б х 0,03) при 1 /гп\ = 0,25: а) динамика доли отказов; б) загруженность 3-го прибора
Входные данные — число обслуживающих приборов, емкость каждого накопителя, а также производительность каждого прибора для случая детерминированной диспетчеризации входящих заявок.
Различимое состояние СМО представлено одномерным массивом 5 = Б(х1, Х-2, • • • , Хк',У1,У2, • • • ) У к) длины 2 к. Здесь х^ = 0,1, где 0 соответствует свободному, а 1 — занятому прибору; уг = 0, т, — наполненность накопителя г-того прибора.
На первом этапе формируем двумерный массив размерности 2А'П'=1 (Уг + 1)) х2к всех возможных комбинаций х% и у г исходя из условий
^Xi = 0,1;
[О ^ Уг ^ ГЩ.
На втором этапе для выделения допустимых состояний 5 из массива всех состояний СМО удалим строки, не отвечающие условию
О
О, т,-
= 0,
Хг = 1,
так как у свободного прибора очереди нет, после чего получим массив размерности п х 2к, где п — число допустимых состояний СМО. При этом каждой строке массива ставим в соответствие строковую переменную, содержащую обозначение данного состояния в нотификации [1, табл. 1]. При детерминированной диспетчеризации в таблице допустимых состояний вводится дополнительный столбец с рассчитанным параметром
Г, = Е
г=1
№
На следующем этапе для каждого номера ] € 1, п найдём те из оставшихся п — 1 состояний, в которые система перейдёт из состояния Sj при поступлении очередной заявки (прямая дуга орграфа состояний) или при выбытии заявки из системы (обратная дуга орграфа состояний).
Рассмотрим два произвольных различных допустимых состояния СМО и соответствующие им вершины орграфа в1, в" € в. Обозначим через в" элементы соответствующих одномерных массивов. Тогда Б1 связано направленной дугой с Б" тогда и только тогда, когда существует такой единственный номер что ф в", \з'1 — = 1. При этом, если з'1 — = 1, то дуга пря-
мая, если же з'1 — з\ = — 1, то дуга обратная. Указанное условие в алгоритме реализовано равенством (в" — = 1. При этом для заполнения матрицы
смежности знак разности в" — не важен.
Матрица смежности для СМО с диспетчеризацией формируется сначала без учета диспетчеризации (1), а затем удаляются лишние элементы выше главной диагонали, соответствующие прямым дугам орграфа состояний системы. Таким образом, для каждого 5' = 5*, \/^ > 5* существует один и только один элемент матрицы смежности Маф-(,§,£) = 1, для которого
Т&1 = пип Т*.
Для оптимального расположения вершин орграфа состояний СМО (рис. 3) рассчитаем их координаты по следующему алгоритму:
у- = Л"4+1;
а* = 2> + ^ Ъ+1 = \о^х^х<+1.
%— 1 V,
Отметим, что данный алгоритм визуализации орграфа оптимален лишь для исследуемых СМО с раздельными накопителями и неоднородными приборами.
5.(100000) 6.(101000) 7.(110000) 8.(111000)
Рис. 3. Орграф СМО сигнатуры Г = (0 х 0,24; 0 х 0,5; 0 х 0,26) с диспетчеризацией
Заключение. Разработанный программный комплекс позволяет получать численные характеристики и автоматически строить граф состояний для практически распространенного случая СМО с различимыми каналами, т. е. СМО с приборами различной производительности и/или раздельными накопителями. При этом нет ограничений на пуассоновское распределение интервалов времени между событиями потоков системы, а также отсутствие последействия, что важно при наличии регулярной составляющей у потоков входящих и обработанных заявок.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. А. 77. Котенко, М. Б. Букаренко, “Аналитическое описание систем массового обслуживания с использованием колец вычетов” / В сб.: Труды седьмой Всероссийской научной конференции с международным участием. Часть 2 / Математическое моделирование и краевые задачи. Самара: СамГТУ, 2010. С. 136-140. [А. P. Kotenko, М. В. Bukarenko, “Analytical description of queuing systems with use residue rings” / In: Proceedings of the Seventh All-Russian Scientific Conference with international participation. Part 2 / Matem. Mod. Kraev. Zadachi. Samara: Samara State Technical Univ., 2010. Pp. 136-140].
2. A. 77. Котенко, М. Б. Букаренко, “Система массового обслуживания с различимыми каналами как конечный автомат” // Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2012. №3(28). С. 114-124. [А. P. Kotenko, М. В. Bukarenko, “A queuing system with distinct devices as the finite state machine” // Vestn. Samar. Cos. Tekhn. Univ. Ser. Fiz.-Mat. Nauki, 2012. no. 3(28). Pp. 114-124].
3. A. 77. Котенко, М. Б. Букаренко, “Анализ высоковолатильных рынков с использованием метода Берга и фильтров Чебышева II рода и статистическое моделирование риска убыточности его инструментов”// Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2011. №3(24). С. 189-192. [А P. Kotenko, М. В. Bukarenko, “High volatile markets analysis with using Berg method and Chebyshev type II filters, and statistical modeling of the risk of loss for its tools” // Vestn. Samar. Cos. Tekhn. Univ. Ser. Fiz.-Mat. Nauki, 2011. no. 3(24). Pp. 189-192].
4. М. Б. Букаренко, “Статистическое моделирование входящего потока требований системы массового обслуживания, включающего детерминированную компоненту” / В сб.: Труды восьмой Всероссийской научной конференции с международным участием. Часть 2 / Математическое моделирование и краевые задачи. Самара: СамГТУ, 2011. С. 136-137. [М. В. Bukarenko, “Statistical modeling of the incoming flow of requests queuing system which includes a deterministic component” / In: Proceedings of the Eighth All-Russian Scientific Conference with international participation. Part 2 / Matem. Mod. Kraev. Zadachi. Samara: Samara State Technical Univ., 2011. Pp. 136-137].
Поступила в редакцию 02/III/2013; в окончательном варианте — 05/V/2013.
MSC: 60К25, 62Р20; 90В22, 68М20, 62М20
SOFTWARE SYSTEM FOR SIMULATION OF STUDY OF QUEUING SYSTEM WITH HETEROGENEOUS SERVERS AND INDEPENDENT QUEUES
A. P. Kotenko, М. B. Bukarenko
Samara State Technical University,
244, Molodogvardeyskaya St., Samara, 443100, Russia.
E-mails: akol959@mail. ru, maxim. bukarenko@gmail. com
Software package for simulation of study and automated visualization of the state or-graph of queuing system with heterogeneous servers and independent limited queues is described. The base of the algorithm is the method of description of such queues with the use of Mealy finite-state machines.
Key words: queuing system, non-homogeneous servers, separate queues, graph visualization, incident matrix of the state orgraph.
Original article submitted 02/111/2013; revision submitted 05/V/2013.
Andrey P. Kotenko (Ph.D. (Phys. & Math.)), Associate Professor, Dept, of Applied Mathematics & Computer Science. Maksim B. Bukarenko, Postgraduate Student, Dept, of of Applied Mathematics & Computer Science.
