Имитационное моделирование экономических процессов

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

a) не гарантируется типичность данных;

b) длительность моделируемого процесса ограничена длительностью реального;

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

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

Физическим или программным методом формируется случайное число Ui, равномерно распределенное на интервале [0,1).

Выполняется программный переход от Ui к случайному числу Xi, имеющему требуемое распределение.

Равномерно распределенное на [0,1) случайное число представляется в компьютере в виде n - разрядной последовательности нулей и единиц, причем в каждом разряде нуль или единица должны наблюдаться с вероятностью 0,5.

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

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

Числа, выдаваемые программными генераторами, называются “псевдослучайными”. В работе Лемера (1951г) формулируется идея использования псевдослучайных чисел следующим образом: ”последовательность псевдослучайных чисел можно считать случайной, если ее каждый элемент непредсказуем для непосвященных и она удовлетворяет ряду традиционных статистических тестов”.

В n - разрядной сетке компьютера количество различных двоичных чисел равно 2n, т.е. ограничено. Поэтому последовательность генерируемых псевдослучайных чисел рано или поздно начнет повторяться. Количество чисел, генерируемых до начала повторения последовательности, называется периодом генератора L. В лучшем случае период генератора L=2n. Тогда говорят, что генератор имеет полный период. Однако для многих генераторов период гораздо меньше, чем 2n , и может зависеть от начального значения генератора.

Хороший программный генератор случайных чисел должен обладать следующими свойствами:

1. Получаемые числа должны быть равномерно распределены в интервале [0,1) и не должны иметь корреляции друг с другом. Данное свойство проверяется различными методами: например, исследуется равномерность заполнения r - мерного единичного гиперкуба точками (U1,U2,…,Ur),(Ur+1,Ur+2,…,U2r),…Существуют и другие тесты проверки равномерности и отсутствия корреляции.

2. Генератор должен обладать быстродействием и не требовать больших затрат памяти (удобство использования на практике).

3. В генераторе должен быть предусмотрен простой способ получения отдельных потоков случайных чисел, т.е. независимых последовательностей чисел. Каждый поток можно выделить для конкретного источника случайности при моделировании. Например, один поток - для интервала прихода транзактов, а другой поток - для интервалов обслуживания. Либо один поток - для интервалов прихода транзактов одного типа, а второй - для транзактов другого типа. В этом случае обеспечивается сравнимость результатов моделирования.

4. Генератор должен иметь достаточно большую длину периода.

4.2 Линейные конгруентные генераторы

Сегодня очень часто применяются линейные конгруентные генераторы, созданные Лемером. В них последовательность целых чисел Z1,Z2,… определяется по формуле

Zi=(aZi-1+c)(mod m)

Где m - модуль (деление по модулю означает взятие остатка от деления на m), a - множитель, с - приращение и Z0 - начальное значение являются неотрицательными целыми числами.

Очевидно, что остаток от деления всегда меньше делителя, поэтому . Чтобы получить случайные числа в интервале [0,1), примем Ui=Zi/m.

Если с=0, то линейный конгруентный генератор называется мультипликативным, а при c>0 - смешанным. Мультипликативные генераторы появились раньше и использовались более интенсивно. Большинство генераторов, которые применяются сегодня, являются мультипликативными. У них не может быть полного периода, однако они более изучены и факт улучшения эффективности в связи с введением смешанных генераторов пока не доказан.

Пусть n - число битов в слове компьютера, которые действительно доступны для хранения данных. Например, если в компьютере используются 32 - битные слова, а крайний левый бит является знаковым то n=31. Модуль стараются выбрать равным m=2n=231. В этом случае избегают явного деление на m, которое заменяется обработкой ситуации переполнения разрядной сетки.

В качестве параметра a обычно выбирают числа вида 8k+3 или 8k+5. В этом случае для мультипликативного генератора можно получить период 2n-2, причем операция умножения заменяется сдвигом и сложением (умножение на степень двойки - это, по сути, сдвиг числа в разрядной сетке).

Например, в одном из предлагаемых в литературе генераторов a=1073741821, с=0, m=231, Z0=65539. Используются 32-разрядные числа, причем первый разряд интерпретируется как знаковый. Полный период такого генератора равен 231-2.

Потоки случайных чисел (независимые последовательности) в линейном конгруентном генераторе задаются посредством определения начального числа каждого потока. Например, если нам нужны потоки в 100000 чисел каждый, то зададим для первого потока Z0 =произвольное значение, для второго потока возьмем начальное значение Z100000 , а для третьего - начальное значение Z200000 и т.д. Т.е., в качестве потоков выступают непересекающиеся смежные последовательности случайных чисел.

В системе GPSS World используется мультипликативный конгруентный генератор с периодом 231-2. Он выдает равномерно распределенные числа в интервале [1, 2 147 483 647], которые затем нормируются и преобразуются в числа из интервала 0 - 0,999999 включительно либо в числа от 0 до 999, в зависимости от контекста использования.

Можно использовать любое количество генераторов без специального их объявления, обращаясь к стандартному числовому атрибуту (системной переменной) RNj. Т.е., RNj - это значение генератора случайных чисел, равномерно распределенных на интервале [0,1). Номер генератора j является одновременно и его начальным значением (Z0). Таким образом, каждый генератор задает отдельный поток случайных чисел. Изменить начальное значение генератора в процессе его работы можно только для первых 7 генераторов (с номерами 1 ч 7) с помощью оператора

RMULT A,B,C,D,E,F,G

Где A, B, …, G - начальные значения генераторов 1, 2, …, 7 соответственно.

Например, RMULT ,,111

Генератор номер 3 получает начальное значение 111. (остальные значения не изменяются).

4.3 Генерирование непрерывных случайных величин согласно заданному закону распределения

Существует много методов генерирования случайных величин. Для всех из них требуется источник независимых и равномерно распределенных в интервале [0,1) случайных чисел Ui (записывается это распределение U(0,1)). Числа Ui программно преобразуются в числа Xi, имеющие заданный закон распределения (например, экспоненциальное распределение, нормальное и т.д.). Мы рассмотрим один из методов такого преобразования - метод обратной функции.

Предположим, что необходимо генерировать случайную величину X, являющуюся непрерывной и имеющую функцию распределения F(x) непрерывную и строго возрастающую, причем для любого x 0<F(x)<1. (Напомним, что функция распределения в точке x равна вероятности события, что случайное число X примет значение меньшее или равное x. Т.е. F(x)=P(X?x)). Например, функция распределения может иметь вид, показанный на рисунке. Обозначим F-1 - это обратная функция к F. Тогда алгоритм генерирования случайной величины X будет следующий:

Генерируется Ui из U(0,1).

Возвращается Xi=F-1(Ui)

Очевидно, что F-1(Ui) всегда будет определено, т.к. 0?Ui<1, а областью значений функции F(x) является отрезок [0,1]. На рисунке 15 этот алгоритм изображен графически.

Рис.15. Графическая иллюстрация метода обратной функции

В системе GPSS World обратное преобразование числа, которое дает генератор RNj, в число с заданным законом распределения выполняется с помощью функций языка PLUS. Язык PLUS является расширением языка GPSS. Эти функции являются стандартными в системе GPSS World и могут быть использованы без всяких предварительных описаний.

Имеются 24 вероятностных распределения. Рассмотрим следующие:

Экспоненциальное (показательное) распределение:

Real=Exponential (RNj, m, s)

Здесь m - смещение распределения, а s - масштабный параметр.

Среднее v=m+s, а дисперсия D=s2.

Нормальное (гауссово) распределение:

Real=NORMAL (RNj, m, s)

Здесь среднее v=m и D=s2.

Параметр RNj этих функций есть номер генератора равномерных случайных чисел.

Тема 5. Анализ выходных данных

5.1 Переходный период стохастического процесса

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

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

Пусть величины Y1, Y2, … представляют стохастический (случайный) процесс, развивающийся во времени. Например, Yj - это средняя длина очереди в течении j - го часа (либо коэффициент использования устройства в течение j-х суток и т.п.) Продолжительность одного прогона на модели равна m. Тогда в результате одного i-го прогона получаем реализацию этих случайных величин yi1, yi2, …, yim. Допустим, что выполняется n прогонов, каждый с разной последовательностью случайных чисел. Т.е. предполагается, что перед началом каждого прогона статистические счетчики переводились в исходное состояние, использовались одни и те же начальные условия и генераторы случайных чисел выдавали разные последовательности величин. Тогда получаем матрицу реализаций {yij} (i=1,…n; j=1,…,m).

y11, …, y1i, …, y1m

y21, …, y2i, …, y2m

……………………..

yn1, …., yni, …, ynm

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

.

Поведение стохастического процесса можно разделить на переходный и установившийся (стационарный) период. Поведение в переходный период зависит от заданных начальных условий I. Например, если моделируется работа универсального магазина, то процесс изменения длины очереди в кассу зависит от числа покупателей, которые изначально находятся в торговом зале. Длительность самого переходного периода также зависит от начальных условий. Функции плотности распределения для каждого столбца матрицы {yij} являются различными.

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

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

В случае, когда ставится задача исследования стационарного режима, возможны три подхода:

А) Разработать модель так, чтобы условия функционирования были типичными с самого начала. Этот подход требует от разработчика знания типичных условий работы и умения внести в модель эти условия, что не всегда реально.

Б) Проводить моделирование так долго, что данные, собранные в начале, “утонут” среди типичных данных и не окажут заметного влияния. Этот подход нежелателен, поскольку долгое моделирование требует значительных финансовых затрат.

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

Момент входа системы в стационарный режим можно определить по графику некоторой основной ее характеристики (обычно используется коэффициент загрузки устройства или памяти). Например, на рис.16 показан график изменения коэффициента загрузки многоканального устройства, из которого очевидно, что стационарный период наступает приблизительно в 1500 единиц модельного времени.

Рис.16. График коэффициента нагрузки многоканального устройства

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

5.2 Оценка средних значений показателей

Допустим, необходимо получить оценку установившегося среднего v=E(Y) процесса Y1, Y2, …. Существуют шесть основных подходов к решению этой проблемы. Рассмотрим один из них - метод репликации и удаления.

Для получения оценок используются наблюдения, не входящие в переходный период работы l в ходе каждого прогона имитационной модели. Допустим, мы выполнили n повторных прогонов имитационной модели, длина каждого прогона равна m наблюдений, причем m существенно превышает длину переходного периода l. Обозначим yij - значение исследуемого показателя, которое было получено в ходе i-го независимого прогона в момент j (как и ранее). Определим величины Xi - среднее по времени i - го прогона, которое соответствует стационарному режиму:

Эти величины являются независимыми и одинаково распределенными. Для них можно найти точечную оценку математического ожидания по формуле:

(5.1)

А также оценку дисперсии:

5.2)

Доверительный интервал, в который попадает математическое ожидание величины Xi c вероятностью 1-б определяется по формуле:

(5.3)

Где - критическое значение распределения Стьюдента с n-1 степенью свободы и вероятностью выхода за интервал б . Может быть вычислено по таблице либо с помощью функции Excel СТЬЮДРАСПОБР(). Вероятность попадания случайной величины в доверительный интервал P=1-б называется надежностью интервальной оценки.

5.3 Получение заданной точности и расчет числа повторных прогонов

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

Выполняется предварительно n прогонов модели и по формулам (5.1) - (5.3) рассчитываются точечные оценки и доверительный интервал при заданном уровне надежности б.

Если половина доверительного интервала, полученного по формуле (5.3), больше в, то рассчитывается приблизительное требуемое число прогонов по формуле

(5.4)

Выполняется дополнительно прогонов и снова рассчитываются характеристики по формулам (5.1) - (5.3).

Размещено на Allbest.ru