Software system for simulation of study of queuing system with heterogeneous servers and independent queues

Full Text

Введение. В нотификации [1] рассмотрим систему массового обслуживания (СМО) с k > 0 неоднородными приборами и раздельными очередями сигнатуры T = T (m1 × µ1 ; m2 × µ2 ; . . . ; mk × µk ), где µi — пропускная способность, mi – число мест в очереди i-того прибора, i ∈ 1, k. Её граф уже не будет графом процесса гибели–размножения, но будет определяться протоколом диспетчеризации входных заявок между приборами. В распространенных на практике случаях аналитические методы моделирования данных СМО не работают. В первую очередь, это связано с тем, что они распространяются только на системы с простейшими потоками событий. Уравнения Колмогорова для марковских процессов составляются при условии существования предела pij (t) , ∆t→0 ∆t λij = lim где pij — вероятность перехода системы из i-того в j-тое состояние. Это влечёт за собой пуассоновское распределение интервалов времени между событиями потока. Кроме того, моделирование СМО с использованием тем или иным образом уравнений Колмогорова подразумевает отсутствие последействия, когда для любых двух непересекающихся отрезков времени число событий, попавших на один из них, зависит только от длины отрезка, но не зависит от числа событий, попавших на другой отрезок. Функциональная связь между моментами появления новых заявок в СМО с самым простым регулятором нарушает данное условие. 50 Комплекс программ имитационного моделирования работы СМО. . . При этом уравнения Колмогорова эффективно решаются только в случае однородных марковских цепей, когда гарантированно существует стационарное решение. Наконец, немаловажна громоздкость системы уравнений Колмогорова для СМО с различимыми каналами, порождённая большим числом состояний системы. Представление такой СМО конечным автоматом (КА) Мили позволило создать комплекс программ для моделирования работы системы с произвольным распределением интервалов времени между появлением входящих/обработанных заявок, который строит матрицу смежности и граф состояний данной СМО как без диспетчеризации (т. е. распределение входящих заявок по доступным приборам произвольно и не зависит от их загруженности), так и в соответствии со следующим протоколом диспетчеризации. Пусть очередная заявка обнаруживает систему в состоянии (x1 , x2 , . . . , xk ; y1 , y2 , . . . , yk ), не являющемся состоянием отказа (1, 1, . . . ,1; m1 , m2 , . . . , mk ). Если существует единственный прибор (с номером i), способный принять заявку (0 yi |xi mi − 1|), то заявка направляется к нему. В противном случае выбираем i-тый прибор, способный обработать заявку с минимальным средним суммарным временем T обслуживания напрвленных к нему заявок: T = min i:0 yi |xi mi −1| µ−1 (yi + 1 + χi ) , i (1) где случайная величина χi ∈ [0; 1] характеризует незавершённость обработки заявки i-тым прибором в момент поступления новой заявки. Программный комплекс содержит следующие программы 1–3. 1. Программа имитационного моделирования потоков входящих и обработанных заявок. В основе работы программы лежит алгоритм [3, 4]. При имитации реальной СМО искусственные реализации входящего и выходящего потоков необходимого объема традиционно предполагают лишь пуассоновское распределение интервалов времени между входящими и обработанными заявками. С помощью КА можно сохранить не только произвольные законы распределения, но и структуру входного/выходного сигнала, под которым понимается эмпирический числовой ряд [T1 , T2 , . . . , Ti , . . . , Tk ] из интервалов времени Ti между i-той и (i − 1)-й по счёту буквой входного алфавита (1 — для входящего потока, или 0 — для потока обработанных заявок) [2]. Применив к этому ряду методы цифровой обработки сигналов, получим оценку необходимого числа искусственных реализаций. На первом этапе выделим характерные циклы моделируемых потоков входящих/обработанных заявок: например, для абонентской активности оператора сотовой связи периоды циклов будут составлять день, неделю, год и т.п. Для этого методом максимальной энтропии Берга оценим спектральную плотность мощности (СПМ) каждого эмпирического сигнала. Выбор порядка p авторегрессионной модели производится по критерию длин минимального описания [3]: M LD[p] = N ln ρp + p ln N, где ai — параметры модели Берга, N — длина сигнала, ρi = ρi−1 (1 − |ai [p]|2 ) монотонно уменьшаются с ростом порядка модели. 51 А. П. К о т е н к о, М. Б. Б у к а р е н к о Выбор частоты отсечки производится по впадинам графика СПМ так, чтобы в полосу пропускания попали все доминирующие циклы. Далее сигнал очищается медианной фильтрацией от импульсных шумов, когда перепады значений сигнала велики по сравнению с дисперсией аддитивного белого шума. Выходной временной ряд yk скользящего медианного фильтра с окном ширины 2n + 1 для k-того отсчета формируется из входного эмпирического ряда (. . . , xk−1 , xk , xk+1 , . . .) формулой yk := med (xk−n , xk−n+1 , . . . , xk−1 , xk , xk+1 , . . . , xk+n−1 , xk+n ) , где med (x1 , . . . , x2n+1 ) := xn+1 при ранжировании значений попавшего в окно вариационного ряда: x(1) = min (x1 , x2 , . . . , x2n+1 ) x(2) ... x(2n+1) = max (x1 , x2 , . . . , x2n+1 ) . Если эмпирический сигнал является суммой Ti = yi + ξi детерминированной и стохастической величин, то сохранение его структуры при имитационном моделировании означает представление каждой искусственной реализации суммой Ti = yi + ξi , где ξi — случайная величина с тем же спектром дисперсий, что и ξi . При этом в качестве тренда yi может выступать как отклик низкочастотного линейного фильтра (например, в соответствии с описанной выше процедурой медианной фильтрации), так и заданная аналитически функция (найденная, к примеру, методами асинхронного гармонического анализа). Моделирование стохастической составляющей ξ(t) производится по спектральному разложению ∞ Uk cos ξ(t) = k=0 kπt kπt , + Vk sin T T при некоррелированности центрированных случайных величин Uk , Vk , дисперсии которых для одинаковых индексов k составляют величины D0 = 1 T T kξ (τ )dτ, Dk = 0 2 T T kξ (τ ) cos 0 kπτ dτ, T k 1, где kξ — автокорреляционная функция ξ(t) на (−T ; T ): ∞ kξ (τ ) = Dk cos k=0 kπτ . T Данный метод позволяет моделировать работу СМО при наличии неслучайных составляющих в интенсивностях входного/выходного потоков, что нарушает условие отсутствия последействия. Программа создает требуемое число искусственных реализаций входящего потока заявок и потоков обработанных заявок для каждого отдельного прибора СМО. 2. Программа имитационного моделирования работы СМО. Эта часть программного комплекса реализована на языке Matlab следующим образом. 52 Комплекс программ имитационного моделирования работы СМО. . . Период моделирования разбивается на оптимальное с точки зрения соблюдения условия ординарности и удобства анализа результатов число тактов. В каждый такт вырабатывается n + 1 сигнал, где n — число приборов системы. Состояние системы определим вектором S1×n = {x1 + y1 , x2 + y2 , . . . , xi + yi , . . . , xn + yn } ; 0 yi mi ; xi ∈ [0, 1], в соответствии с принятой нотификацией. Каждый сигнал прихода/обработки заявки является двоичным числом sij . Описание сигналов представлено в таблице. Описание сигналов имитационной модели Порядковый номер сигнала в такте Значение sij j=1 1 j=1 0 j = 2, n + 1 1 j = 2, n + 1 0 Описание Появление входящей заявки. Если St = m1 + 1, m2 + 1, . . ., mi + 1, . . ., mn + 1, то St+1 = St , заявка получает отказ. В противном случае она направляется на свободный прибор в соответствии с принятой диспетчеризацией. St+1 = St Выработка j-тым прибором обработанной заявки: xj+1 + yj+1 = xj + yj − 1, xj + yj 0, xj + yj = 0 1, St+1 = St На основе разработанной нотификации и представления СМО конечным автоматом Мили на языке Matlab была разработана система имитационного моделирования работы СМО с различимыми каналами для произвольных параметров потоков входящих и обработанных заявок. Она позволяет учитывать следующие характеристики (рис. 1): распределение интервалов времени между появлениями входящих заявок; интенсивность потока входящих заявок (для пуассоновского распределения) и математическое ожидание интервала времени между входящими заявками (для прочих типов распределения); дисперсию потока входящих заявок; число приборов системы (произвольное конечное); распределение интервалов времени между обработанными заявками, покидающими систему; интенсивности потоков обработанных заявок (для пуассоновского распределения) и математические ожидания интервала времени между обработанными заявками (для прочих типов распределения) для каждого прибора системы; дисперсии потоков обработанных заявок для каждого прибора; емкости накопителей приборов (произвольные конечные). Программа позволяет получить набор состояний исследуемой СМО на каждом такте в течение заданного пользователем времени работы имитационной модели. Все время моделирования T необходимо так разбить на N равных тактов времени τ , чтобы выполнялись следующие условия:N τ T; P (n > 1, τ ) = 0. Второе условие означает, что вероятность появления более одной входящей заявки за один такт равна нулю, что обеспечивает ординарность потока входящих заявок. Аналогичным образом соблюдается ординарность и для потоков обработанных заявок. 53 А. П. К о т е н к о, М. Б. Б у к а р е н к о Для выполнения первого условия сначала полагаем длину вектора генерируемых псевдослучайных чисел равной N1 = T /Mλ [z], затем увеличиваем значение N1 на единицу до тех пор, пока сумма сгенерированных псевдослучайных чисел не окажется больше времени работы модели T : N1 xi T . При i этом длина такта τi = k mini zi , где k < 1 и задано пользователем. Отметим, что минимально необходимый для ординарности потока промежуток времени между сигналами выработки заявки τ out = min {τiout } i∈[1, n] Рис. 1. Окно ввода параметров может быть меньше продолжительности такта τ in , определяемого исходя из характеристик потока входящих заявок. В этом случае полагаем τ = min{τ in , τ out }. Разработанная программа имитационного моделирования СМО с различимыми каналами позволяет найти следующие характеристики системы: динамика загруженности j-того прибора; доля отказов; динамика вероятности простоя. Средняя загруженность j-того прибора рассчитывается как отношение математического ожидания числа заявок, обслуженных j-тым прибором за данный промежуток времени, к максимально возможному числу заявок, которое прибор в состоянии обработать за данный временной интервал: (n+1)ki+1 Lcp = (n + 1)−1 (kj+1 − kj )−1 i=(n+1)k1 xij + yij , mj + 1 где kj — номер такта. Доля отказов рассчитывается как отношение числа заявок, заставших систему в состоянии St = {m1 + 1, m2 + 1, . . . , mi + 1, . . . , mn + 1} и пришедших kj+l за время i=kj , l = const, к суммарному числу заявок, пришедших в систему за указанный интервал времени. Пример динамики значений данных характеристик представлен на рис. 2. 3. Автоматическое построение графа состояний СМО. Разработанная программа построения графа состояний СМО решает следующие задачи: построение списка допустимых состояний; расчет и визуализация матрицы смежности вершин; расчет координат вершин для наиболее эффективной визуализации орграфа состояний СМО; визуализация графа. Ее интерфейс реализован на VBA for MS Excel, а визуализация осуществляется специально разработанным приложением Graph Theory Toolbox для Matlab. Обмен данными между MS Excel и Matlab происходит посредством специального набора инструментов Excel Link Toolbox, позволяющего вызывать функции Matlab из модуля VBA for MS Excel. 54 Комплекс программ имитационного моделирования работы СМО. . . б a Рис. 2. Результаты моделирования СМО сигнатуры T = (2 × 0,15; 4 × 0,05; 6 × 0,03) при 1/mλ = 0,25: а) динамика доли отказов; б) загруженность 3-го прибора Входные данные — число обслуживающих приборов, емкость каждого накопителя, а также производительность каждого прибора для случая детерминированной диспетчеризации входящих заявок. Различимое состояние СМО представлено одномерным массивом S = S(x1 , x2 , . . . , xk ; y1 , y2 , . . . , yk ) длины 2k. Здесь xi = 0, 1, где 0 соответствует свободному, а 1 — занятому прибору; yi = 0, mi — наполненность накопителя i-того прибора. На первом этапе формируем двумерный массив размерности k k 2 i=1 (yi + 1) ×2k всех возможных комбинаций xi и yi исходя из условий xi = 0, 1; 0 y i mi . На втором этапе для выделения допустимых состояний S из массива всех состояний СМО удалим строки, не отвечающие условию yi = 0 ← xi = 0, 0, mi ← xi = 1, так как у свободного прибора очереди нет, после чего получим массив размерности n × 2k, где n — число допустимых состояний СМО. При этом каждой строке массива ставим в соответствие строковую переменную, содержащую обозначение данного состояния в нотификации [1, табл. 1]. При детерминированной диспетчеризации в таблице допустимых состояний вводится дополнительный столбец с рассчитанным параметром k Tj = i=1 xij + yij + 2χij , µi j = 1, n. На следующем этапе для каждого номера j ∈ 1, n найдём те из оставшихся n−1 состояний, в которые система перейдёт из состояния Sj при поступлении очередной заявки (прямая дуга орграфа состояний) или при выбытии заявки из системы (обратная дуга орграфа состояний). 55 А. П. К о т е н к о, М. Б. Б у к а р е н к о Рассмотрим два произвольных различных допустимых состояния СМО и соответствующие им вершины орграфа S , S ∈ S. Обозначим через si , si элементы соответствующих одномерных массивов. Тогда S связано направленной дугой с S тогда и только тогда, когда существует такой единственный номер i, что si = si , |si − si | = 1. При этом, если si − si = 1, то дуга прямая, если же si − si = −1, то дуга обратная. Указанное условие в алгоритме реализовано равенством (si − si )2 = 1. При этом для заполнения матрицы смежности знак разности si − si не важен. Матрица смежности для СМО с диспетчеризацией формируется сначала без учета диспетчеризации (1), а затем удаляются лишние элементы выше главной диагонали, соответствующие прямым дугам орграфа состояний системы. Таким образом, для каждого s = S, ∀t > S существует один и только один элемент матрицы смежности Madj (s, t) = 1, для которого Tst = min Tt . ∀t>s=S Для оптимального расположения вершин орграфа состояний СМО (рис. 3) рассчитаем их координаты по следующему алгоритму: k (xi + yi ); Xt = i=1 Yt+1 = Yt + ← Xt = Xt+1 ; 0 ← Xt = Xt+1 . Отметим, что данный алгоритм визуализации орграфа оптимален лишь для исследуемых СМО с раздельными накопителями и неоднородными приборами. Рис. 3. Орграф СМО сигнатуры T = (0 × 0,24; 0 × 0,5; 0 × 0,26) с диспетчеризацией Заключение. Разработанный программный комплекс позволяет получать численные характеристики и автоматически строить граф состояний для практически распространенного случая СМО с различимыми каналами, т. е. СМО с приборами различной производительности и/или раздельными накопителями. При этом нет ограничений на пуассоновское распределение интервалов времени между событиями потоков системы, а также отсутствие последействия, что важно при наличии регулярной составляющей у потоков входящих и обработанных заявок. 56 Комплекс программ имитационного моделирования работы СМО. . .

About the authors

Andrey P Kotenko

Samara State Technical University

Email: ako1959@mail.ru
(Ph. D. (Phys. & Math.)), Associate Professor, Dept. of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russia

Maksim B Bukarenko

Samara State Technical University

Email: maxim.bukarenko@gmail.com
Postgraduate Student, Dept. of of Applied Mathematics & Computer Science 244, Molodogvardeyskaya st., Samara, 443100, Russia

References

Supplementary files