6. ОПРЕДЕЛЕНИЕ ХАРАКТЕРИСТИК СТАЦИОНАРНОГО СЛУЧАЙНОГО ПРОЦЕССА МЕТОДОМ МОНТЕ-КАРЛО ПО ОДНОЙ РЕАЛИЗАЦИИ

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

Рис. 8.22

В качестве примера рассмотрим работу -канальной системы массового обслуживания с местами в очереди, граф состояний которой показан на рис. 8.22. Предположим, что поток заявок, переводящий систему из состояния в состояние (слева направо), — стационарный, но не пуассоновский, например поток Пальма с произвольным законом распределения интервала времени Т между заявками. Время обслуживания одной заявки тоже распределено не по показательному закону, а по произвольному закону . Так как процесс, протекающий в системе, немарковский (потоки событий — непуассоновские), то его не удается описать с помощью стандартного математического аппарата марковских случайных процессов — обыкновенных дифференциальных уравнений для вероятностей состояний и алгебраических уравнений — для предельных вероятностей состояний. И вообще, попытка описать данный случайный процесс с помощью аналитических зависимостей привела бы к чересчур громоздкому аппарату, не оправдывающему себя на практике. Единственным практически пригодным методом исследования подобных немарковских систем является моделирование процесса методом Монте-Карло.

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

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

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

Рис. 8.23

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

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

Если же из состояния система перешла не в то она будет циркулировать не по состояниям а по состояниям вероятности этих состояний тоже будут стремиться к постоянным:

но уже к другим, чем

Таким образом, в приведенном примере процесс, протекающий в системе, будет стационарным, но не эргодическим, и его вероятностные характеристики существенно зависят от начального периода (начального поведения системы). Ясно, что моделирование такого процесса с помощью одной (хотя бы и очень длинной) реализации недостаточно для получения его вероятностных характеристик.

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

Если система имеет бесконечное множество возможных состояний, то, мы знаем, даже при стационарности всех потоков событий, предельного режима при может не существовать — мы это видели примере системы массового обслуживания с неограниченной очередью (см. § 6 гл. 5), где при очередь при растет неограниченно. Однако, если предельный режим существует, то при моделировании процесса методом Монте-Карло можно ограничиться одной реализацией.

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

Для пояснения изложенного рассмотрим пример, относящийся к моделированию методом Монте-Карло работы немарковской системы массового обслуживания с очередью.

Пример. Имеется двухканальная СМО с очередью. Число мест в очереди заявка, пришедшая в момент, когда все три места в очереди заняты, получает отказ и покидает систему. Поток заявок — пальмовский, т. е. интервалы времени между заявками представляют собой независимые случайные величины, распределенные по одному и тому же (непоказательному) закону (рис. 8.24). Время обслуживания одной заявки — также случайная величина, распределенная по непоказательному закону (рис. 8.25), отличному от но одинаковому для всех заявок.

Требуется, моделируя работу СМО методом Монте-Карло и располагая только одной длинной реализацией, оценить приближенно предельные характеристики системы (при ):

— вероятности состояний (вероятности того, что будут заняты 0,1, 2 канала; вероятности того, что в очереди будут находиться 0, 1, 2, 3 заявки);

— среднее число занятых каналов;

— среднее время ожидания заявки в очереди; дисперсию времени ожидания заявки в очереди;

— вероятность отказа (того, что заявка покинет СМО необслуженной).

Рис. 8.24

Рис. 8.25

Построить схему моделирования и схему обработки его результатов.

Решение. Граф состояний системы имеет вид, показанный на рис. 8.26. Число состояний конечно; из каждого состояния можно перейти в каждое; потоки событий, переводящие систему из состояния в состояние, стационарны (хотя и непуассоновские); из этого заключаем, что система обладает эргодическим свойством и моделирование по одной реализации возможно.

Рис. 8.26

Предположим для простоты, что в начальный момент система находится в состоянии (свободна). Начнем моделирование с того, что разыграем на оси поток заявок, т. е. ряд случайных точек — моментов прихода соответствующих заявок — первой, второй и т. п. (рис. 8.27). Розыгрыш потока заявок производится следующим образом. Строится функция распределения случайной величины Т — интервала между заявками:

и разыгрывается значение случайной величины как это описано в § 2; для этого берется функция, обратная F, от случайного числа R от 0 до

Расстояние откладывается от начала координат; получается момент прихода первой заявки. Затем процедура розыгрыша повторяется (разумеется, уже при другом R) и новое значение откладывается от получается момент прихода второй заявки, и т. д. Таким образом мы построим цепочку моментов прихода заявок (рис. 8.27). Разумеется, эту цепочку надо сделать достаточно длинной, разыграв, во всяком случае, не менее нескольких сот значений случайной величины Т.

Рис. 8.27

Рис. 8.28

Изобразим процедуру моделирования с помощью наглядной схемы (рис. 8.28). Вверху мы поместим ось времени (0) с отмеченными на ней моментами поступления заявок. Ниже ее мы поместим еще пять осей: (1), (2), (3), (4), (5). На осях (1) и (2) мы будем изображать состояния первого и второго каналов (жирная черта — «занят», тонкая — «свободен»). На осях (3), (4), (5) мы будем изображать состояния первого, второго и третьего мест в очереди (жирная черта — «занято», тонкая — «свободно»). Все пять осей имеют тот же отсчет времени, что и ось (0).

До момента — прихода первой заявки — все каналы и все места в очереди свободны.

В момент приходит первая заявка и занимает первый канал. Сколько времени он будет занят — решается розыгрышем. Для этого мы подвергнем случайное число R (разумеется, новое) преобразованию , где — функция распределения времени обслуживания:

Первое разыгранное значение времени обслуживания обозначаем и откладываем на оси (1) от точки с абсциссой отмечая его жирной линией (рис. 8.28). В момент прихода второй заявки первый канал еще занят; заявка занимает второй канал. Разыграем еще одно значение , обозначим его и отложим жирной линией на оси (2) от точки с абсциссой

Заявка , пришедшая в момент, когда оба канала заняты, становится в очередь, занимает в ней первое место (ось ) и ждет до того момента, когда освободится один из каналов. В нашем случае раньше освобождается канал (2) — в этот момент точка с оси (3) перескакивает на ось (2) — и снова разыгрывается время обслуживания этой заявки. На оси (2) строится новый жирный участок, а ось (3) продолжается тонкой линией — место в очереди свободно.

Не будем продолжать подробное описание процедуры розыгрыша реализации — она достаточно ясна из рис. 8.28. На этом рисунке против каждого участка занятости канала (места в очереди) для удобства обработки проставлен номер заявки, занимающей это место; можно проследить, как заявка «путешествует» с последних мест в очереди на первые, затем — на обслуживание. Заявка, получившая отказ, отмечается звездочкой (она покидает СМО необслуженной).

Предположим что моделирование реализации продолжено нами достаточно долго (настолько долго, что влияние начальных условий уже перестает сказываться). Посмотрим, как по этой реализации определить интересующие нас вероятностные характеристики работы СМО. Вероятности того, что будут заняты 0, 1, 2 канала, найдем следующим образом. Разделим всю ось на участки соответственно числу занятых каналов. Участки времени, на которых не занят ни один канал, отметим цифрой 0, один канал — цифрой 1, два канала — цифрой 2. На большом участке времени Т сложим длины всех участков, помеченных нулем — получим сумма длин всех участков, помеченных единицей, будет двойкой —

Очевидно,

При большом Т вероятности будут приближенно равны отношениям соответствующих времен к общему времени:

Заметим, что участок Т целесообразно отсчитывать не с самого начала процесса, где еще сказывается влияние начальных условий, а от более удаленного по времени момента 0, где влияние начальных условий уже практически перестает сказываться.

Найдем вероятности того, что в очереди будут стоять 0, 1, 2, 3 заявки. Снова разобьем большой участок оси времени Т на части, помеченные 0, 1, 2, 3, на которых в очереди находится соответственно 0, 1, 2, 3 заявки. Складывая длины всех одинаково помеченных участков и деля суммы на Т, получим:

Среднее число занятых каналов получится обычным способом как математическое ожидание дискретной случайной величины Z — числа занятых каналов:

Среднее время ожидания заявки в очереди находим следующим образом: рассмотрим ряд заявок, поступивших на большом участке времени Т в моменты

и для каждой из них непосредственно подсчитаем время ожидания в очереди равное нулю, если заявка была сразу принята к обслуживанию (или получила отказ), и сумме времен ожидания этой заявки на разных осях ((3), (4) и (5)), если она стояла в очереди. Среднее время ожидания заявки в очереди приближенно найдется как среднее арифметическое этих времен:

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

Дисперсия времени ожидания найдется аналогичным образом, как среднее арифметическое квадратов времен ожидания минус квадрат среднего времени ожидания:

Наконец, вероятность отказа найдется на большом участке времени Т как отношение числа N заявок, помеченных звездочкой (получивших отказ), к общему числу N заявок, поступивших за это время: