Введение Метод моделирования различных сложных систем, прежде всего с использованием средств вычислительной техники, занимает важное место при анализе различных реальных объектов [1, 2]. Применительно к социотехническим системам, основное содержание которых определяется процессом взаимодействия человека и технических систем, этот метод приобретает особую актуальность, поскольку поведение одного из ключевых компонентов социотехнических систем - человека - обычно не поддается точному описанию, для описания его поведения используются различные стохастические, неопределенные категории и понятия [3, 4]. Аналитическое описание подобных элементов эффективно при использовании их стохастического описания, т. е. на основе аппарата теории вероятностей. Часто одним из основных направлений взаимодействия является попытка различных субъектов воспользоваться одним и тем же ограниченным ресурсом. Возникает конкуренция за право воспользоваться ресурсом. Для разрешения возникшего противоречия используются различные аналитические и эвристические методы, в частности системы массового обслуживания. Точный анализ подобных систем зачастую практически невозможен, поэтому для их изучения используются методы моделирования. Многие из подобных систем массового обслуживания могут быть эффективно описаны на основе более узкого их класса - многоканальных систем массового обслуживания. Таким образом, одним из важных элементов моделей социотехнических систем являются многоканальные системы массового обслуживания, поэтому эффективность процедуры моделирования искомых систем существенно зависит от эффективности моделирования требуемых многоканальных систем массового обслуживания, а также от точности воспроизведения ими реальных объектов. Основная их особенность - индивидуальный, не совпадающий с другими схожими субъектами, характер всех характеристик каждого участвующего в работе социотехнической системы субъекта, что вызвано разнородностью их места и роли в составе системы. При этом весь поток заявок порождается конечным числом источников, и каждая порождаемая заявка имеет свои индивидуальные характеристики. Существующие процедуры моделирования многоканальных систем массового обслуживания [5] не в полной мере позволяют учесть индивидуальные особенности реальных социотехнических систем. Также эти процедуры недостаточно приемлемы с точки зрения скорости работы программ, реализующих их. В работе приводится процедура моделирования многоканальных систем, которая позволяет в определенной мере устранить указанные недостатки. Для повышения быстродействия предлагается заменять исходные функции распределения их приближениями, которые могут быть достаточно быстро реализованы. Среди близких работ укажем на работу [6]. Таким образом, одной из основных особенностей предлагаемой ниже процедуры моделирования многоканальной системы массового обслуживания является индивидуализированность процесса моделирования до уровня отдельного вызова и отдельного обсуживающего прибора. То есть все поступающие вызовы и каждый обслуживающий прибор могут иметь свои индивидуальные характеристики (например, функции распределения длительностей обслуживания, интервалов поступлений и т. п.), и предлагаемая процедура позволяет учесть эти особенности. Заметим, что подавляющее большинство существующих систем моделирования обычно объединяют вызовы и приборы в классы однотипных, и в процессе моделирования выбор любого из них не приводит к отличиям в процессе вычислений. Рекуррентные соотношения Пусть имеется система массового обслуживания (СМО), состоящая из некоторого количества k обслуживающих приборов. Ограничений на длину очереди нет. Предположим, что одновременное окончание обслуживания на нескольких приборах невозможно, т. е. имеет вероятность нуль. Данное условие, в частности, выполняется в случае, когда функции распределения длительностей обслуживания непрерывны. Введем обозначения: - - промежуток времени между поступлениями в систему (n -1)-го и n-го вызовов; - - длительность обслуживания n-го по порядку поступления вызова; - tn - момент n-го освобождения одного из обсуживающих приборов; - pn - длина очереди в момент tn; - - номера вызовов в порядке поступления их в систему (т. е. при i < j), которые в момент tn находятся в очереди при условии, что pn > 0; при pn = 0 множество Q(n) не определено; - jn - номер обслуживающего прибора, который освобождается в момент tn (1 ≤ jn ≤ r); - () - время дообслуживания вызова на приборе с номером i в момент времени tn + 0, т. е. с учетом того, что при наличии очереди на jn-й прибор поступит из очереди вызов в соответствии с некоторой политикой выбора вызовов из очереди; - ρ(n) и ( при i < j) - число и номера свободных приборов в момент времени tn; если все приборы заняты, то ρ(n) = 0 и I(n) = Ø. Тогда . Необходимо также описать процедуры выбора вызовов из очереди и выбора прибора для обслуживания из числа свободных приборов в ситуациях, когда свободны более одного прибора. Первая из этих процедур задается с помощью специальной функции. Именно, задана функция переключения , определяющая номер вызова из находящихся в очереди, который выбирается для обслуживания в момент времени tn. Отметим, что для классических дисциплин выбора FIFO (выбор в порядке поступления), LIFO (выбор последнего из пришедших из числа находящихся в очереди) и RS (случайный выбор из очереди) функция переключения соответственно равна , и , где Rand(N) - это функция случайного выбора числа из набора {1, 2, …, N}. Отметим, что в этом случае в момент tn на jn-й прибор поступает вызов с номером , и . Далее, справедливо соотношение (1) где - номер первого после момента tn поступающего в систему вызова, поскольку в момент tn СМО уже покинуло n вызовов, в очереди находится pn вызовов, и на обслуживании находятся r - ρ(n) вызовов; δ(n) - номер последнего до момента tn+1 поступившего в систему вызова. При этом необходимо заново перенумеровать все вызовы в Qn+1: (2) Отметим, что последний случай может отсутствовать, если pn > pn+1; при этом справедливо равенство Нетрудно убедиться в справедливости равенства (3) Для выбора обсуживающего прибора в момент tn в случае, если имеется несколько свободных приборов, задана функция выбора прибора . В простейшем случае можно взять (т. е. выбрать первый свободный) либо (т. е. обслуживающий прибор выбирается случайным образом из числа свободных). Тогда можно записать соотношение При этом, как и выше было выполнено для Qn+1, необходимо перенумеровать номера приборов в I(n +1). Процедура формирования I(n+ 1) описана в следующем разделе. Составим рекуррентные соотношения для вычисления всех введенных величин. Справедливы равенства (4) Здесь мы воспользовались тем, что если , т. е. нет очереди, то все поступившие до момента tn вызовы обслужены; но т. к. число освобождений обслуживаний в момент tn равно n, то, следовательно, следующим по номеру поступающим вызовом, который и займет свободный прибор, является (n + 1)-й вызов. При этом в момент tn+1 либо освобождается один из уже занятых в момент tn приборов (в этом случае в системе освободится вызов через время ), либо один из поступающих вызовов занимает свободный прибор и освобождается раньше остальных - для i-го поступающего вызова: их номера изменяются от (n + 1) до (n + ρ(n)) эта величина равна . Процедура вычислений Полученные выше рекуррентные соотношения для всех основных характеристик СМО в моменты окончаний обслуживаний достаточно сложно взаимосвязаны между собой. В частности, в вычислении почти всех характеристик в момент tn необходимо знание величины tn+1; в свою очередь, величина tn+1 находится на основе значений многих характеристик в момент tn. Поэтому необходимо описать последовательность вычисления всех характеристик на основе полученных рекуррентных соотношений, которая бы позволила обойти указанный круг взаимных зависимостей. Как показано ниже, полученные ранее в работе соотношения позволяют решить эту задачу. То есть они позволяют реализовать следующую рекуррентную процедуру вычисления всех текущих характеристик системы, перечисленных выше, в момент tn+1, при известных значениях тех же характеристик в момент tn и значениях входных величин { n ≥ 1}. Именно, пусть для некоторого n известны значения величин p(n), Qn, I(n), (). Тогда последовательно вычисляются следующие величины, с опорой на значения известных или уже вычисленных величин. 1. Находим jn: если в момент tn система пуста, то полагаем jn = 0, в противном случае полагаем (см. (4)). 2. Если pn > 0 (т. е. в момент tn в очереди есть вызовы и, следовательно, I(n+1)=Ø), то, очевидно, jn > 0, и полагаем , - промежуток времени, через который произойдет следующее завершение обслуживания в системе, . Затем вычисляем , (см. (3)) и (см. (1)). В множестве I(n+ 1) производится перенумерация вызовов на основе соотношения (2). Также полагаем , . Процедура завершается. 3. Если pn = 0 (т. е. в момент tn в очереди нет вызовов и, следовательно, возможно I(n+1) = Ø), то для вычисления I(n+1), () выполняется следующая циклическая процедура. Полагаем k = 1, I0(n+1) = I(n)\ jn, S0 = 0, () и переходим к шагу 4. Перенумеруем элементы множества I0(n + 1) следующим образом: пусть и I0(n + 1) = {}. Тогда 4. Для данного значения k полагаем Sk = Sk-1 + . Пусть k < d(0), т. е. еще имеются свободные приборы. Полагаем , для полагаем , и . При k ≥ d(0) (свободных приборов нет) имеем Ik(n+ 1) = Ø, и аналогично случаю pn > 0 имеем (). В случае Ik-1(n+1) ≠ Ø перенумеруем элементы множества Ik(n+1) следующим образом: пусть Ik-1(n+1) = {}, причем αi,k-1 < αj,k-1 при i < j. Тогда d(k) = d(k-1) - 1 и 5. Если k < d(0), то проверяется выполнение условий: Sk-1 < σ(k), Sk ≥ σ(k). Если данные условия выполнены, то в этом случае процесс нахождения момента tn+1 закончен: полагаем tn+1 = tn + σ(k-1), (),= k-1; (поступивший после момента tn k-й по порядку вызов занял прибор с номером ). Далее, (см. (1)). На этом процедура вычисления характеристик системы в момент tn+1 завершается. Если условия не выполнены, т. е. Sk < σ(k), то увеличиваем k на единицу (т. е. k: = k+1) и переходим к шагу 2. Процедура вычисления характеристик СМО в момент tn+1 закончена. Как следует из содержания приведенной процедуры, на каждом шаге используются либо заданные (исходные) данные, либо уже вычисленные до данного текущего момента вычислений. Аппроксимация исходных данных Полученные выше выражения и приведенная рекуррентная процедура позволяют сделать ряд важных выводов. Прежде всего, они позволяют находить последовательность основных характеристик СМО в произвольный момент времени, причем исходные случайные характеристики могут иметь любые распределения и даже быть зависимыми. Далее, описанная в работе СМО непрерывна в схеме серий с точки зрения сходимости по вероятности [4, 5]. То есть если исходные характеристики системы {; ; n ≥ 1} сходятся по вероятности к некоторому предельному набору {; ; n ≥ 1}, то все текущие характеристики системы, включая длины очередей {pn; n ≥ 1}, также сходятся к соответствующим значениям [8-10]. Это следует из того, что в описанной выше процедуре использованы функции min(), x(), argmin() на конечном множестве значений, условные операторы типа «если …, то …», которые сводятся к индикаторной функции x() (см., например, (2)), и непрерывны по отношению сходимости по вероятности [5-7]. Отсюда вытекает возможность использования различных аппроксимаций исходных данных с целью упрощения процесса вычислений, в частности при имитационном моделировании рассмотренных в работе моделей. Рассмотрим данный вопрос более подробно. Наибольший практический интерес представляет случай, когда все случайные величины {; ; n ≥ 1} независимы. Более того, с точки зрения практических применений можно исходные распределения указанных случайных величин заменить на их приближенные выражения, характеристики которых можно было бы достаточно быстро вычислять. Примерами подобных приближений являются, например, ступенчатые функции распределения, соответствующие приближениям исходных распределений с помощью дискретных распределений. Особо важное их достоинство - это возможность использования эмпирических функций распределения, построенных на основе статистических данных, - эти функции как раз и имеют ступенчатый вид. Отметим, что любую функцию распределения можно приблизить ступенчатыми функциями с любой заданной степенью точности и заданной сеткой дискретных значений функции. Еще одним классом функций, представляющим интерес с рассматриваемой точки зрения, являются смеси показательных распределений, поскольку показательные распределения обладают рядом важных достоинств: показательное распределение обладает свойством отсутствия памяти, когда будущее поведение случайной величины не зависит от ее прошлого, что часто упрощает анализ систем; характеристики показательного распределения вычисляются на основе достаточно простых выражений. Однако класс показательных распределений не замкнут, поскольку смеси показательных распределений представляют собой гамма-распределения. Поэтому целесообразно рассмотреть более широкий класс распределений - класс гамма-распре-делений, который замкнут, т. е. сумма и предел любой последовательности гамма-распре-делений, если он существует и не вырожден, также является гамма-распределением. Возникает вопрос о полноте класса гамма-распределений и смесей показательных распределений, т. е. о возможности приближения любого распределения с помощью распределений из указанных классов с любой заданной точностью, и о процедуре получения непосредственно требуемых приближенных выражений для заданных функций распределения. Также необходимо описать процедуру нахождения искомого приближения для заданной функции распределения или заданного набора статистических данных. Приведем формализованную постановку задачи. Дана произвольная функция распределения G(x) неотрицательной случайной величины, т. е. G(x) = 0 при x < 0, и заданная малая величина ε требуемой точности приближения. Необходимо: 1) показать, что существует набор показательных функций распределения и вероятностей , - таких, что для любого x > 0 выполнено (5) 2) найти параметры и , - такие, при которых выполнено (5). Решение указанной задачи рассмотрено автором в работе [11], где приведен один из вариантов требуемого алгоритма. Пусть - требуемая точность приближения. 1. Разобьем интервал [0; 1] возможных значений функции G(x) на промежутки длиной меньше ε/2. Для этого достаточно выбрать и разбить [0; 1] на промежутки для , . Находим числа , k < N. Отметим, что если G(x) - непрерывная функция, то . В результате ось аргументов Ox разобьется на промежутки . 2. Для каждого фиксированного k, находим минимальное значение M = M(k) - такое, что где - функция гамма-распределения с параметрами (m, α), Г(α) - гамма-функция Эйлера, Из следствия 4 [11] вытекает, что последнее неравенство при достаточно больших M выполняется. Искомое значение M может быть найдено методами прямого перебора, деления отрезка пополам, другими методами направленного перебора. Вариант процедуры нахождения M(k) приведен в разделе 3 работы [11]. 3. В качестве результата берется выражение Нетрудно показать, что для всех x. Использование при моделировании распределений Gε(x) вместо функций распределений G(x) позволит значимо повысить скорость формирования искомых случайных величин, поскольку имеются достаточно эффективные алгоритмы моделирования гамма-распределений. Заключение В работе для многоканальной системы массового обслуживания, в которой все вызовы имеют индивидуальные характеристики поступления и обслуживания в соответствии с их особенностями в искомой социотехнической системе и переключения с обслуживания одного вызова на другой осуществляются в соответствии с заданной функцией переключения, определяемой принятой в социотехнической системе политикой разрешения конфликтных ситуаций, получены рекуррентные соотношения для длин очередей, перечня свободных приборов, перечня номеров вызовов, ожидающих обслуживания и ряда других характеристик в последовательные моменты окончаний обслуживаний вызовов. Описана процедура последовательного вычисления всех указанных характеристик системы с учетом их достаточно сложной взаимосвязи. Для повышения эффективности процесса практической реализации процедуры моделирования в работе приведена процедура приближения исходных распределений с помощью смесей гамма-распределений.