3.2. Метод Монте-Карло

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

Идея метода Монте-Карло состоит в том, что вместо аналитического описания СМО производится “розыгрыш” случайного процесса, проходящего в СМО, с помощью специально организованной процедуры. В результате такого получается каждый раз новая, отличная от других реализация случайного процесса. Это множество реализаций можно использовать как некий искусственно полученный статистический материал, который обрабатывается обычными методами математической статистики. После такой обработки могут быть получены приближенно любые характеристики обслуживания.

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

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

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

Метод Монте-Карло (метод статистических испытаний) включает следующие этапы:

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

2. Установление распределения вероятностей для стохастических переменных.

3. Установление интервала случайных чисел для каждой стохастической переменной и генерация случайных чисел.

4. Имитация поведения системы путем проведения многих испытаний и получение оценки моделируемой характеристики системы при фиксированных значениях параметров управления. Оценка точности результата.

Описание этапов:

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

где F, R - вектор-функции;

I_i, U_i, S_i - векторы соответственно входных, выходных переменных и переменных состояния системы в тактовый момент моделирования i.

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

Существует два подхода в имитационном моделировании динамических процессов.

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

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

Второй этап. Случайные величины, используемые в ИМ, могут быть дискретными или непрерывными. В первом случае необходимо знать их распределения, во втором - плотности распределений. Эти зависимости могут быть известны из теории, определены в результате специальных исследований либо заданы в качестве гипотезы. Точность модели (при прочих равных условиях) зависит от того, насколько точно заданы указанные распределения (плотности распределений).

Третий этап. Моделирование случайных величин при компьютерных имитационных экспериментах производится с помощью датчика псевдослучайных чисел, предусмотренного в любом современном языке программирования. Обычно это датчик случайных чисел с равномерным распределением на интервале [0, 1]. Если известны вероятности наступления событий, то, используя такой датчик, можно отвечать на вопросы: «Какое из n возможных событий произошло?» или «Какое значение приняла случайная величина?»

Предположим, что в ИМ используется случайная величина Х, принимающая дискретные значения х₁, х₂, ..., х_n, с вероятностями соответственно р₁, р₂, ...,р_n, ( = 1). Получение некоторой реализации этой переменной в модели производится следующим образом.

Строится функция распределения случайной величины Х. Указанная функция определяется посредством равенства F(Х)=Σp_k, в котором суммирование распространяется на все индексы, для которых х_k < Х. С помощью датчика случайных чисел получают случайное число u из отрезка [0, 1].

Из равномерности распределения получаемых случайных чисел следует, что вероятность получения случайного числа из произвольного интервала, включенного в [0, 1], равна длине этого интервала. Поэтому вероятность реализации Х=x_k равна вероятности попадания полученного от датчика случайного числа u в произвольный интервал длиной p_k на отрезке [0, 1]. Можно, таким образом, утверждать, что если очередное число u датчика удовлетворяет неравенствам 0<u≤р₁, то имеет место реализация Х = х₁; в случае р₁ < u ≤ p₁ + р₂ - реализация Х= х₂ и т.д. В общем случае для k = 2,..., n: если то Х=х_k.

Границы указанных неравенств совпадают со значениями построенной выше функции распределения F(Х).

Удобнее, однако, иметь дело не с дробными значениями границ интервалов, в которые попадает случайное число u, а с их целочисленными значениями, тем более, что с помощью датчиков случайных чисел можно генерировать числа из любого диапазона. Чтобы получить целые значения границ интервалов, достаточно умножить все р_k на 10^d, где d - целое, минимальное значение которого равно максимальной точности (максимальному числу знаков после десятичной точки) чисел р_k, k = 1, ..., n. Например, если {р_k} = {0,3; 0,153; 0,5; 0,047}, то минимальное значение d равно 3 (все р_k нужно умножить на 1000). Таким образом, 10^d опредёляет длину интервала значений рассматриваемой случайной величины в ИМ.

Четвертый этап. Точность статистических оценок параметров реальной системы зависит от числа наблюдений (объёма выборки). Погрешности в оценках обусловлены как статистическим характером самой модели, так и влиянием начальных данных (начального состояния имитационной системы), а также возможной автокорреляцией последовательных значений некоторого параметра в процессе моделирования. Очевидно, что с увеличением числа испытаний точность моделирования должна возрастать. Ввиду того что увеличение объема выборки связано с ростом затрат на моделирование, важно уметь определять минимальное число испытаний, необходимое для достижения заданной точности оценки с заданной вероятностью.

Широкое распространение получили два метода статистических испытаний. Один из них предполагает проведение достаточно большого числа N последовательных наблюдений в течение одного прогона модели (одного сеанса имитирования).

Другой метод заключается в реализации m независимых прогонов модели, т.е. в m-кратном повторении одного и того же цикла имитирования. При этом, если мы хотим получить в сумме N наблюдений, в течение каждого прогона можно делать по N/m (допустим, что это число целое) наблюдений. Оба метода дают примерно одинаковый результат.

Пусть значения у_t (t = 1, ..., N) представляют собой результаты N последовательных измерений значений случайной величины у во время одного и того же сеанса имитации. Среднее по времени значение у определяется выражением

Обозначим через µ математическое ожидание случайной величины у. Тогда для достаточно большого N получаем

Оценка дисперсии (если временной ряд не является автокоррелированным) имеет вид

где D(у) - дисперсия случайной величины у.

Для оценки качества результатов, полученных методом Монте-Карло при неизвестной дисперсии наблюдаемой случайной величины, предположим, что Z - характеристика, которая должна быть определена (вероятность события, математическое ожидание, дисперсия и т.п.), а ξ - ее значение, уточняемое по мере накопления данных, остающееся случайным вследствие ограниченности числа N проведенных наблюдений. В этих условиях можно говорить о вероятности р(|Z – ξ| < ε) по отношению к интересующей нас характеристике. Величина |Z – ξ| представляет собой погрешность в оценке Z, а ε - некоторый допустимый ее предел.

Из неравенства Чебышёва следует

р(|Z – ξ| < ε) ≥ 1 - D_ξ(N)/ε².

Из этого неравенства следует

D_ξ(N) < (1 - р) ε²,

откуда при заданных р и ε и при известной зависимости D_ξ(N) можно найти предельно необходимое N.

Известно, что истинная дисперсия выборочного распределения для расчетного среднего обратно пропорциональна суммарному числу наблюдений N, т.е.

D_ξ(N) = d/N,

где d не зависит от N.

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

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

Предельное число наблюдений N₀ определяется формулой

N₀ = d ⁄ [(1 – p) ε²].

При любом числе наблюдений больше N₀ обеспечивается требуемая точность.

Пример имитационной модели очереди на разгрузку.

Наблюдения за транспортными потоками в течение 200 дней показали, что число транспортных средств, ежедневно прибывающих на предприятие с грузами колеблется от 0 до 5. Вероятность прихода 0, 1, …, 5 транспортных средств показана в таблице:

Число прибывающих транспортных средств	Частота появления транспорта (количество дней)	Вероятность	Интегральная вероятность (значение функции распределения)	Интервал случайных чисел
0	26	26/200=0,13	0,13	от 01 до 13
1	34	34/200=0,17	0,30	от 14 до 30
2	30	30/200=0,15	0,45	от 31 до 45
3	50	50/200=0,25	0,70	от 46 до 70
4	40	40/200=0,20	0,90	от 71 до 90
5	20	20/200=0,10	1,00	от 91 до 00
Итого	200	1,00	-	-

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

Ежедневный темп разгрузки	Частота разгрузки транспорта (количество дней)	Вероятность	Интегральная вероятность (значение функции распределения)	Интервал случайных чисел
1	10	10/200=0,05	0,05	от 01 до 05
2	30	30/200=0,15	0,20	от 06 до 20
3	100	100/200=0,50	0,70	от 21 до 70
4	40	40/200=0,20	0,90	от 71 до 90
5	20	20/200=0,10	1,00	от 91 до 00
Итого	200	1,00	-	-

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

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

День	Количество транспорта, простаивающее с предыдущего дня	Случайное число RND	Количество прибывшего транспорта за день	Количество транспорта, ожидающего разгрузки	Случайное число RND	Количество разгруженного транспорта
1	-	50	3	3	10	2
2	1	48	3	4	24	3
3	1	61	3	4	03	1
4	3	18	1	4	32	3
5	1	85	4	5	23	3
6	2	23	1	3	59	3
7	-	08	-	-	95	-
8	-	54	3	3	34	3
9	-	17	1	1	34	1
10	-	12	-	-	51	3
11	-	80	4	4	08	2
12	2	69	3	5	48	3
13	2	24	1	3	66	3
14	-	84	4	4	97	4
15	-	92	5	5	03	1
Итого	12 (общий простой)	-	36 (всего прибыло)	-	-	32 (всего разгружено)
Среднее число	12/15=0,8	-	36/15=2,4	-	-	32/15=2,133

Имитационная модель управления запасами. В фирме, реализующей продукцию в течение 300 дней регистрировался дневной спрос. Распределение вероятностей величины спроса указано в таблице:

Спрос на продукцию	Частота спроса (количество дней продаж)	Вероятность	Интегральная вероятность (значение функции распределения)	Интервал случайных чисел
0	15	15/300=0,05	0,05	от 01 до 05
1	30	30/300=0,10	0,15	от 06 до 15
2	60	60/300=0,20	0,35	от 16 до 35
3	120	120/300=0,40	0,75	от 36 до 75
4	45	45/300=0,15	0,90	от 76 до 90
5	30	30/300=0,10	1,00	от 91 до 00
Итого	300	1,00	-	-

Когда менеджеры делают заказ, чтобы возобновить запасы продукции, его выполнение происходит с лагом (задержкой) в 1, 2 или 3 дня. Это означает, что время восстановления запаса подчиняется вероятностному распределению. В следующей таблице указаны сроки, вероятности сроков выполнения заказов и интервалы случайных чисел, которые удалось определить на основе информации о 50 заказах:

Срок выполнения заказа	Частота поставок с данным лагом (количество дней)	Вероятность	Интегральная вероятность (значение функции распределения)	Интервал случайных чисел
1	10	10/50=0,20	0,20	от 01 до 20
2	25	25/50=0,50	0,70	от 21 до 70
3	15	15/50=0,30	1,00	от 71 до 00
Итого	50 заказов	1,00	-	-

Условия резервирования, которые собирается имитировать менеджмент фирмы – делать заказ в объёме 10 единиц продукции при запасе на складе 5 штук (точка восстановления запаса). Каждый заказ обходится в 10 ден.ед. в день, а хранение каждой единицы продукции – в 5 ден.ед. в день, одна упущенная продажа – в 80 ден.ед. Цель эксперимента – оценить величину средних ежедневных затрат для этой системы управления запасами.

Решение. Реализуется четырёхшаговый процесс имитации:

1. Каждый имитируемый день начинается с проверки, поступил ли сделанный заказ. Если заказ выполнен, то текущий запас увеличивается на величину заказа (в данном случае – на 10 единиц).

2. Путём выбора случайного числа генерируется дневной спрос для соответствующего распределения вероятностей.

3. Рассчитывается итоговый запас, равный исходному запасу за вычетом величины спроса. Если запас недостаточен для удовлетворения дневного спроса, спрос удовлетворяется, насколько это возможно. Фиксируется число нереализованных продаж.

4. Определяется, снизился ли запас до точки восстановления, принятой равной 5 единицам. Если да, причем не ожидается поступления заказа, сделанного ранее, то делается заказ.

Результаты имитационного эксперимента представлены в табл. 15.

Ежедневные затраты на заказы = Затраты на один заказ  Среднее число заказов в день = 10 · 0,3 = 3 ден.ед.

Ежедневные затраты на хранение = Затраты на хранение одной единицы в течение дня  Средняя величина конечного запаса = 5 · 3,8 = 19 ден.ед.

Ежедневные упущенные продажи = Прибыль от упущенной продажи  Среднее число упущенных продаж в день = 80 · 0,5 = 40 ден.ед.

Общие ежедневные затраты = 3 + 19 + 40 = 62 ден.ед.

Таблица 15

Результаты имитационного модельного эксперимента

День	Поступление заказа	Начальный запас	RND	Спрос	Конечный запас	Потери продаж	Делать заказ?	RND	Срок выполнения заказа
1	-	10	07	1	9	-	0 - Нет	-	-
2	-	9	60	3	6	-	0 - Нет	-	-
3	-	6	77	4	2	-	1 - Да	52	2
4	-	2	49	3	0	1	0 - Нет	-	-
5	-	0	76	4	0	4	0 - Нет	-	-
6	10	10	95	5	5	-	1 - Да	37	2
7	-	5	51	3	2	-	0 - Нет	-	-
8	-	2	16	2	0	-	0 - Нет	-	-
9	10	10	14	1	9	-	0 - Нет	-	-
10	-	9	85	4	5	-	1 - Да	82	3
Всего		63		30	38	5	3
Среднее значение (дневное)		6,3		3	3,8	0,5	0,3

159