Имитационное моделирование

                                          ,

запишем правые части уравнений  в матрицу D, здесь − элементы матрицы интенсивностей. Интервал интегрирования (0; 10), число разбиений 1000.

(На приведенных графиках вероятности состояний обозначены как yi).

 

Правые части уравнений:

Вектор начальных условий

 

 

Построим график переходного процесса в системе

 

       

4. Моделирование процесса, протекающего в данной системе.

Введем переменный массив sj, элементы которого – суммарное время пребывания системы в данном состоянии j, матрицу В – индикатор состояний (каждый столбец соответствует одному состоянию). Например, при выборе столбца 3 система  находится в состоянии 2.  Напомним, что в этой главе элементы массивов нумеруются так 0,1,2,3…

 

Моделирующая программа.

 

 

 

 Рассмотрим операторы программы по порядку. Задается начальное значение модельного времени t. Вводится матрица  iw соответствующая начальному состоянию системы по индикатору состояний и строится цикл while до достижения времени моделирования tm. Определяется номер состояния k, в котором система находится в текущий момент времени. В цикле вычисляем все времена , через которые система может перейти в другое состояние.                 Находим минимальное из этих времен . Так как цепь Марковская, то

 она удовлетворяет условию  отсутствия последействия, и случайные времена между переходами распределены по показательному закону. Они могут быть найдены с помощью оператора . Здесь 1 показывает, что вычисляется одно значение, а − интенсивность соответствующего перехода. Полученное время суммируется с временем, которое система провела в текущем состоянии s. Определяется номер ind состояния, в которое переходит система, и этот номер присваивается индикатору В. Отношения времени пребывания в каждом состоянии к полному времени моделирования, принимаются за оценки стационарных вероятностей состояний. Для рассмотренного примера и времени моделирования получим

 

 

 

 

 Сравнивая результаты моделирования при различных прогонах с различными числами шагов и точные значения стационарных вероятностей состояний, делаем вывод о хорошей сходимости результатов.

 

5. Моделирование систем массового обслуживания

 

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

 

 5.1 Пошаговое моделирование СМО

 

   Пусть при помощи стандартной функции MathCad runif(n,a,b) получена последовательность из десяти случайных чисел, равномерно распределенных на интервале (0,1).

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

 Время обслуживания показательное  с параметром  μ и 

имеет тот же закон распределения Пусть интенсивность потока , а среднее время обслуживания tob=0.2 часа. Так как в момент t=0 система простаивает, то возможным  событием может быть лишь приход заявки, занятие ей канала и начало обслуживания − событие А. Первая заявка поступит в систему через часа, или в 0.067 часа.

Теперь возможны следующие случаи − в систему поступит очередная  заявка и встанет в очередь (событие  А1) или канал освободится  и система перейдет в состояние простоя − событие В. Какое из них произойдет, определяется временем прихода требования и временем обслуживания. Так как , то заявка поступит в , а , то конец работы канала в . Таким образом, и произойдет событие А1 − вторая пришедшая заявка встанет в очередь. Далее возможны два события А1 или В1 − заявка покинет канал и требование из очереди попадет в канал. Определим время прихода ,   и произойдет событие В1. К числу обслуженных заявок нужно прибавить единицу, вычислить время заявки в очереди . Используем четвертое случайное число для определения времени обслуживания , а время прихода , , происходит событие В, число обслуженных заявок , канал освободится и система будет простаивать до . В этот момент заявка займет канал и покинет его в . Очередное требование поступит в (система в это время простаивает), причем обслуживание продлится до 1.23часа, а в 1.2 часа  произойдет событие А1 − следующее требование  встанет в очередь. В (событие В1) заявка из очереди попадает на обслуживание и покинет канал в 1.54 (по случайному числу ,а пришедшая заявка в1.28 займет место в очереди ( ). Массив случайных чисел исчерпан. Всю последовательность событий и связанных с ними изменениями состоянии системы оформим в виде протокола моделирования.                                        

                                 

                                     Протокол моделирования.

Модельное время	Событие	Число обслуженных заявок	Общее время в очереди	Число заявок в очереди	Состояние системы
0.0	-	0	0	0	Простой
0.067	А	0	0	0	Работа
0.223	А1	0	0	1	Работа
0.244	В1	1	0.021	0	Работа
0.374	В	2	0,021	0	Простой
0.440	А	2	0.021	0	Работа
0.99	В	3	0.021	0	Простой
1.157	А	3	0.021	0	Работа
1.2	A1	3	0.021	1	Работа
1.23	В1	4	0.051	0	Работа
1.28	A1	4	0.256	1	Работа

   

    Легко заметить, что время моделирования (машинное время) никак не связано с «реальным» временем, протекающим в системе. Так первая заявка поступит через 0.067 часа после начала процесса, а машинное время на совершение простейшего расчета этой величины имеет порядок  долей секунды.

       Рассмотрение последовательности событий в приведенном примере дает представление о работе моделирующей программы. Аналогичные построения приведены, в частности в [3,4].

 

 

 

5.2 Моделирование СМО с ограничением на длину очереди.

 

Для дальнейших сравнений результатов, сначала проведем расчет параметров простейшей СМО.    Тогда, - интенсивность входного потока, - среднее время обслуживания и . Число каналов n , каналы нумеруются индексом k . Вероятности  состояний  в стационарном случае находим для дальнейшего сравнения по обычным формулам [4,5], так как моделирование позволяет получить только стационарные значения вероятностей.

 

   

    

 

 

 

 

 

 

Построение моделирующей программы.

     Пусть N – число пришедших заявок, No – число обслуженных, а Np – число потерянных  заявок. (В случае неограниченной очереди Np =  0),

Tob – суммарное время затраченное на обслуживание и Toch – суммарное время пребывания заявок в очереди. Для сокращения записи при работе программы и выводе результатов введем матрицу-столбец rem, тогда

В программе и подпрограммах  для этой матрицы в качестве формальных параметров используются обозначения  re, re1, re2 и т.д. Зададим матрицу для времен постановки заявки в очередь Och (формальный параметр O, O1 …).

Основным алгоритмом моделирования  выберем, так называемый, «просмотр активностей». Этот алгоритм заключается в том, что на каждом шаге моделирования просматриваются времена «механизмов» составляющих систему, и выбирается  тот, время которого на текущий момент модельного времени наименьшее и который  должен произвести действие. В нашем примере моделирования СМО, такими «механизмами» являются: 1) входной поток, 2) каналы обслуживания (по их числу) и 3) очередь. Действиями являются: для входного потока – приход заявки, для канала – начало обслуживания и окончание обслуживания, для очереди – постановка в очередь и уход из очереди в канал. Например, пусть следующая заявка должна поступить в систему в = 12.42, первый канал освободится в 12.56, второй канал свободен, (тогда его время ), третий обслужит заявку в 12.39, а четвертый в 13.05. Минимальное из этих времен 12.39. Следовательно, будет выполнено действие - третий  канал освободится, к числу обслуженных заявок прибавится единица, суммарное время обслуживания Tob увеличится на интервал времени, в течение которого заявка была в канале.       

(Более подробное описание пошагового  моделирования приведено в начале данной главы).

Времена срабатывания «механизмов» записываются в матрицу so (в качестве формальных параметров подпрограмм s, s1 ит.д.).          

Для просмотра активностей используется функция,

которая определяет номер элемента массива x при выполнении условия , где ”=” логическое равенство. Например, операция

определит минимальный элемент  массива s , выделит этот элемент  и присвоит переменной ind  номер этого элемента. Ниже приведена программа моделирования и ее описание. Заметим, что иногда для большей наглядности и удобства работы с программой ее можно разделить на ряд подпрограмм. Например, выбор канала,  в который попадет следующая пришедшая заявка можно осуществить в виде подпрограммы vk

Тогда, при обращении из основной программы к  vk(kz,B), эта подпрограмма определит наличие свободных каналов и укажет его номер kz. Здесь В индикатор занятости канала  Bk = 0, если k-ый канал свободен Bk = 1 – занят.

      Зададим модельное время   t = 0 и число занятых мест в очереди m = 0 – очереди нет. Построим цикл while , который будет выполняться пока  модельное время, не превысит времени моделирования tm. Обсуждение выбора tm приведено в разделе анализа результатов.

     Цикл 

предназначен для исключения возможного зацикливания программы при простое  СМО. В этом фрагменте времена срабатывания каналов приводятся к времени прихода следующей заявки плюс некоторая малая положительная величине eps, которая должна быть  на 3-4 порядка меньше промежутка между поступающими заявками и не оказывать влияния на процесс моделирования. Этой же цели служит строка

 в завершающей части программы. 

           Блок, связанный с поступлением новой заявки в систему. 

     Если переменная ind равна нулю, то активным «механизмом» является входной поток и ближайшим действием будет приход заявки. Сначала проверим, превышено ли предельное значение занятых мест в очереди mo.

   Если очередь равна или  превышает число мест, то фиксируется  приход заявки и ее потеря.  Если очередь не превысила предельное значение, то запоминаем время прихода предыдущей заявки, увеличиваем число пришедших заявок и определяем время прихода следующей.

 

    С помощью подпрограммы vk  находим свободный канал kz , если их нет, то kz = 0, заявка встает в очередь m = m+1 и запоминается время постановки в очередь Om = rt.

    Если свободный канал есть, то находим время обработки заявки ts, назначаем время конца работы этого «механизма», устанавливаем признак занятости этого канала и суммируем время обслуживания Tob = Tob + ts.

 

           Блок, связанный с обслуживанием заявки в канале. 

    Пусть теперь активен один из каналов обслуживания (ind > 0).