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

Федеральное агентство  железнодорожного транспорта

Российской Федерации

Уральский государственный университет  путей сообщения

__________________________________________________________________

Кафедра высшей и прикладной математики

 

 

 

                                                   

                                                                                         

 

 

 

 

 

 

Курсовой проект по дисциплине: «Имитационное моделирование»

 

                            

 

 

 

 

 

 

 

 

Выполнила: Меньшикова А.А. ПИЭ-319

Проверил:  Скачков П.П.

 

 

 

 

 

 

 

 

 

Екатеринбург

2011

СОДЕРЖАНИЕ

 

 Введение.                                                                                                  3

1. Метод Монте-Карло. Решение детерминированных задач.            

 Моделирование задач имеющих стохастическую природу.               5

   2. Случайные числа.                                                                        6   

3. Вероятностно-статистические аспекты метода Монте-Карло  

  и имитационного моделирования (ИМ).                                                11            

4. ИМ Марковских процессов.                                                                14

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

8. Список литературы                                                                               35

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Введение.

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

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

Компьютерное ИМ предполагает выполнение ряда последовательных действий.

1. Описание реальной системы с выделением структуры, динамического взаимодействия элементов, факторов неопределенности и состояний системы, в которых она может находиться.
2. Создание блоковой схемы объекта с указанием  состояний его элементов и возможных переходов между ними.
3. Построение моделирующей программы на специальном языке ИМ или общем языке программирования.
4. Проигрывание различных возможных ситуаций на модели.
5. Верификация модели и программы на основе анализа полученных результатов и их сравнения с теорией процесса и (или) информацией о функционировании реального объекта.

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

,

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

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

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

     Область применения  ИМ в настоящее время можно  разделить на две основные части.

1. Теоретические задачи математики, математической физики, физики, химии…., в частности
• вычисление многомерных интегралов;
• обращение и псевдообращение матриц;
• вычисление констант;
• решение дифференциальных уравнений в частных производных;
• диффузия в смесях.
2. Практические задачи организации управления

• анализ производственных технологических процессов;
• планирование и прогнозирование инвестиционных проектов;
• задачи социально-психологического характера;
• разработка военной тактики и стратегии.

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

 

 

 

 

 

 

 

 

 

                                                                    

 

 

 

 

 

 

 

 

 

 

 

1. Метод  Монте-Карло (МК).

 

       ИМ можно считать  развитием метода МК, разработанного  в 50-х годах прошлого века. Основная  идея  этого метода состоит  в использовании выборок  для получения оценок искомых характеристик изучаемых объектов. Задача, при этом, формулируется таким образом, чтобы алгоритм решения использовал случайные числа соответствующих законов распределения [3,4]. Достаточно сложно представить себе, как формализовать полностью детерминированную задачу (вычисление определенных интегралов, например) для решения ее с помощью выборок. (Существенное  значение при этом имеют методы  получения последовательностей случайных чисел). Этот вопрос детально рассматривается в следующей главе.

   Найти площадь фигуры  методом Монте-Карло, ограниченной  линиями:

U(t):=4t-t² , x1:=1 и x2:=4

 

Решение:

 

 

 

 

 

 

 

 

нахождение числа точек попавших на фигуру!

 

 n=759 – количество точек, находящихся в выделенной фигуре

P=0.753

Sp=9.036

 

 

 

 

               S=9                   

        

 

 

 

 

 

 

 

 

2. Случайные числа (СЧ)  и методы их получения.

 

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

                         

 математическое ожидание М(X)=1/2 и дисперсия D(X)=1/12.

      Наиболее распространенным методом получения такой последовательности является мультипликативная конгруэнция (в различных модификациях) [1,3]. Два целых числа x  и y называются конгруэнтными по модулю m (m –целое число), если |x–y| = km. Таким образом,  y конгруэнтно x по модулю m, если |x–y| делиться на m без остатка. Например, при x = 12589 и m = 10, y = 9 конгруэнтно x по модулю 10, а при x = 1223 и m = 2, y = 1 конгруэнтно x по модулю 2. Метод состоит в получении последовательности по рекуррентной формуле

,  k = 0,1,2,…

где входные параметры. Такой алгоритм приводит к повторению псевдослучайных чисел начиная с некоторого k. Можно доказать [1], что при

a = 100003, и – девятизначном целом нечетном числе не делящимся на 5, получится не повторяющихся случайных чисел. Очевидно, что при решении задачи необходимо, чтобы полученной последовательности было достаточно для прогона модели.

   Продемонстрируем получение  первых трех псевдослучайных  чисел на примере,

пусть

   Тогда, последовательно имеем

и т.д.

Здесь операция mod(Z,V) определяет остаток от деления Z на V. 

Число неповторяющихся случайных чисел можно существенно увеличить следующим простым способом. После выбора вычисляются значения и вырезаются числа стоящие, например, в 11,12 и 13 разрядах. Пусть в этих разрядах оказались числа 2, 0, 7 тогда 0,207 следующее случайное число.  Этот метод может использоваться и самостоятельно для небольших выборок. После этого находим , описанным выше способом мультипликативной конгруэнции, и получаем последующую серию из чисел. Угол α выбирается произвольно, но меньше 0,5 угловой секунды.

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

     Тест на равномерность последовательности проводится по обычной схеме обработки опытных данных. Интервал (0;1) делится на k частей, определяются частоты и по критерию согласия (например, Пирсона) принимаем гипотезу о равномерном  законе с некоторым уровнем значимости.

       Тест на стохастичность обычно проводят по методу серий. При этом вся исследуемая последовательность делится на элементы первого и второго рода

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

…aaabbbbaabaaabbbbbabab…

Так как случайные числа в  этой последовательности предполагаются независимыми и равномерно распределенными на интервале (0, 1), то теоретическая вероятность появления серии длиной j в последовательности длиной L в N опытах определяется формулой Бернулли

.

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

 

Проверка стохастичности чисел

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

  

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 Проверка независимости элементов последовательности псевдослучайных чисел {xi} проводится на основании вычисления корреляционного момента. Две случайные величины называются независимыми, если закон распределения каждой из них не зависит от того, какое значение приняла другая величина. Введем в рассмотрение случайную величину , где τ − величина сдвига последовательности. Корреляционный момент двух случайных величин X и Y   с реализациями и определяется по формуле

,

 

 

Проверка на независимость.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Ниже приведена сводка результатов  при различных величинах сдвига t.

 

 

Условие независимости  выполняется плохо

 

 

    Общим подходом построения последовательностей случайных чисел с произвольными законами распределения является метод обратной функции. Пусть требуется создать выборку СЧ, имеющих закон распределения . Если случайное число, из последовательности, имеющей равномерное распределение на интервале (0,1), то возможное значение непрерывной случайной величины X c заданной функцией распределения  , является корнем уравнения

  или   .

  Пример 1. Непрерывная случайная величина Х распределена по показательному закону, заданному функцией распределения

.

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

.

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

    Пример 2. Непрерывная случайная  величина Х распределена по  закону равномерной плотности на интервале (a,b), заданному функцией распределения

.

Подставляя функцию в уравнение  и разрешая его относительно , получим

.

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

,

где m и σ математическое ожидание и среднее квадратичное отклонение получаемой  случайной величины, − СЧ равномерно распределенные на интервале (0,1). Для практических целей достаточно просуммировать 12 таких чисел (n=12).   

    

3. Вероятностно-статистические аспекты  метода Монте-Карло и имитационного моделирования

 

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

N=10. Меняя продолжительность прогона от 100 до 10000 точек, вычислим средние выборочные и дисперсии для каждого из этих случаев. Для наглядности приведем программу вычислений площади круга еще раз. Результаты расчетов сведем в таблицу.

 

 

 

 

 

 

 

 

 

 

 

 

 

 

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

1. С ростом числа генерируемых точек (продолжительности прогона модели) оценки площади круга приближаются к точному значению.

 

 

 

 

 

 

 

 

 

 

 

 

   

 

 

 

 

 

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

2. Прогоны модели, отличающиеся друг от друга  только последовательностями случайных чисел, дают различные оценки при одном и том же значении n.
3. Отметим, что влияние переходных условий уменьшается при усреднении результатов прогонов (это хорошо видно из строки таблицы для средних выборочных).
4. Дисперсия рассматриваемой случайной величины существенно уменьшается при увеличении продолжительности прогона с 100 до 200, а затем ее уменьшение незначительно. Этот факт так же характерен для имитационных моделей, что позволяет подобрать оптимальное для рассматриваемей задачи значение n, по критерию точность− затраты машинного времени.

 

 

 

 

 

 

 

 

 

 

 

 

4.Моделирование Марковских процессов.

4.1 Марковская цепь с дискретным временем.

 

     Рассмотрим, для простоты, Марковскую цепь с тремя состояниями (k=3). (Подробное изложение теории Марковских цепей см. [4,5]). Переходы между состояниями происходят мгновенно в фиксированные моменты времени. Вероятности переходов из любого состояния Si в любое другое Sj считаются заданными и равны pij. Моделирование системы требует симуляции стохастических воздействий на систему (случайных чисел и случайных событий ). Получение СЧ описано в главе 3, а случайные события, образующие полную группу, с заданными вероятностями можно моделировать следующим образом. Пусть событие А имеет вероятность р(А), а противоположное событие − (1-р(А)) . Разыграем одну реализацию случайной величины имеющей равномерное распределение на интервале (0,1). Если  , то полагают, что произошло событие А, если , то произошло противоположное событие. Приведем пример с большим числом событий. Пусть события А,В,С образуют полную группу и попарно несовместны, и р(А)=0.3, В – р(В)=0.5 и С– р(С)=0.2. Тогда, если полученное случайное число х  меньше 0.3, то полагают, что произошло событие А , если  х  меньше 0.8, но больше 0.3 – событие В и при х больше 0.8 – событие С. При многократном повторении элементарных  опытов по определению вероятности событий  по приведенной схеме частота появления  каждого события будет стремиться к его вероятности.

   Поставим конкретную задачу. Для Марковской цепи с тремя состояниями задана матрица вероятностей переходов за один шаг  Р.

1. Составить размеченный граф состояний этой Марковской цепи, определить, является ли цепь регулярной.

2. Найти стационарное распределение  вероятностей состояний. 

3. Выполнить моделирование системы  и сравнить полученные результаты с результатами, полученными ранее в пункте 2.

Решение. 1. Составим граф состояний.

 

                 1/3                                              0

                                             1/3

                                        1/2

 

                 1/3                  0

                                                   1/2         1/2

 

 

                                           1/2

          

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

2. По формулам, [5] найдем стационарное распределение вероятностей. Запишем систему алгебраических уравнений соответствующую данному матричному

 

          

         Þ        .              

 

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

Решим систему уравнений:

                    Þ               ,

(0,231; 0,461; 0,308).

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

    Примем, что в начальный момент времени система находится в состоянии S0  и Q(0) = (1, 0, 0). Пусть число шагов моделирования или продолжительность прогона равна ns. Введем матрицу В − индикатор состояния, (например, столбец показывает, что система после последнего шага находится в состоянии S0)  и вектор so для суммирования числа попаданий в каждое из состояний. Так как в начальный момент времени система находилась в состоянии , то , а и .

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

jm− счетчик числа шагов, x −  случайные числа, равномерно распределенные на интервале (0,1). Значение индекса k определяет номер состояния, из которого выходит система, а номер состояния, в которое  осуществляется переход, получим из соотношений  переходных вероятностей. Элементы вектора sо,  есть числа попаданий системы в данное состояние и они будут возрастать на 1, как только система попадает в это состояние. Формальными параметрами программы являются число шагов N и вектор s.

   Краткое описание работы программы.

   Счетчику jm присваивается начальное значение 0, выбирается  столбец соответствующий начальному состоянию и строится цикл while до достижения значения N счетчика  jm. Внутри цикла определяется значение индекса k или номера состояния, из которого будет проходить переход, затем, в соответствии со значениями переходных вероятностей, находится состояние i , в которое перейдет система. В рассматриваемом примере, вероятности переходов из состояния в равны . Тогда, если следующее случайное число х меньше  или равно 1/3, то произойдет событие , если , то и при , и число возрастет на единицу. Индикатор состояний приводится в положение соответствующее совершенному переходу. К счетчику шагов прибавляется единица. Далее расчет продолжается аналогично. Результатом  моделирования будут компоненты вектора s − числа попаданий в каждое из состояний. Точные (вычисленные по формулам [5]) стационарные значения вероятностей состояний представим в виде матрицы

_.

 

Программа, моделирующая процесс, протекающий  в цепи, приведена ниже.

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

                                                            

 

 

 

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

    

 

4.2 Марковская цепь с непрерывным временем

 

Рассмотрим систему с k  состояниями. Переходы между состояниями происходят мгновенно в случайные моменты времени. Вероятности переходов из любого состояния Si в любое другое Sj являются функциями от времени pij(t). Если  случайный процесс, протекающий в системе, обладает свойством отсутствия последействия, то говорят, что задана Марковская цепь с непрерывным временем. Интенсивностью перехода из состояния Si в состояние Sj называется предел ,

где -вероятность перехода на интервале времени .

             Матрица вероятностей состояний удовлетворяет системе дифференциальных уравнений Колмогорова. Эта система может быть записана в матричной

                                                                                   

или в координатной форме

 

                                             

 

Здесь -матрица, составленная из производных , i = 1, 2, 3 вероятностей состояний в момент времени t. Для решения системы (1.2.5) или (1.2.6) необходимо, как обычно, задать начальные условия:

                                ,   ,            

где , i = 1, 2, 3 -заданные числа, причем

Определение 1.2.2. Распределение вероятностей называется стационарным, если вероятности состояний не зависят от времени или .

Для стационарного распределения  вероятностей системы  , тогда из (1.2.5) получаем систему алгебраических уравнений

,

или в координатной форме:

                                                    

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

   Рассмотрим, для примера,  Марковскую цепь с тремя состояниями.  Пусть задана матрица интенсивностей переходов Λ и  начальное распределение вероятностей состояний

 

 

Требуется:

1. Составить размеченный граф состояний этой Марковской цепи, определить, является ли цепь регулярной.

  2. Найти стационарное распределение  вероятностей состояний.

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

4. Выполнить моделирование системы и сравнить полученные результаты моделирования с результатами, полученными в пункте 2.

 

Решение.

 

1.   Составим граф состояний.

 

 

                                                2

                                           3

 

                         3             1             4           4

 

 

                                                                                                   

                

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

2. По формулам найдем стационарное  распределение вероятностей:

                                                 , 

 

Тогда стационарное распределение  вероятностей состояний Sq определим [5]

           3. Решим  уравнения Колмогорова в системе  MathCAD при помощи стандартной функции Rkadapt при начальных  условиях (в начальный момент времени система находится в состоянии S0 ). Сведем систему к двум уравнениям, используя условие нормировки