46362

АНАЛИТИЧЕСКИЕ И ИМИТАЦИОННЫЕ МОДЕЛИ

Книга

Информатика, кибернетика и программирование

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

Русский

2013-11-21

7.07 MB

29 чел.

PAGE  262

МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ

ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ

Технологический институт

Федерального государственного образовательного учреждения высшего профессионального образования

«Южный федеральный университет»

ПРИОРИТЕТНЫЙ НАЦИОНАЛЬНЫЙ ПРОЕКТ «ОБРАЗОВАНИЕ»

В.И. ФИHАЕВ

Е.Н. ПАВЛЕНКО

Е.В. ЗАРГАРЯН

АНАЛИТИЧЕСКИЕ

И ИМИТАЦИОННЫЕ

МОДЕЛИ

Допущено Учебно-методическим объединением вузов по образованию в области автоматизированного машиностроения (УМО АМ) в качестве учебного пособия для студентов высших учебных заведений, обучающихся по направлению подготовки дипломированных специалистов – «Автоматизированные технологии и производства» (специальность 210200 – «Автоматизация технологических процессов и производств (в энергетике))».

Таганрог 2007


УДК 518.5.001.57(075.8)

В.И.Финаев, Е.Н. Павленко, Е.В.Заргарян. Аналитичес-кие и имитационные модели: Учебное пособие.  Таганрог: Изд-во Технологического института ЮФУ, 2007. - 310 с.

ISBN 978-5-8327-0268-1

Учебное поcобие пpедназначено для cтудентов высших учебных заведений, изучающиx дисциплины «Моделиpование cиcтем», «Компьютерное моделирование». В учебном пособии приведены основные теоретические положения и методы имитационного моделирования.

Табл. 8. Ил. 30. Библиогр.: 22 назв.

Рецензенты

Целых А.Н., д-р. техн. наук, профессор, директор регионального (областного) центра новых информационных технологий, проректор по информатике ТТИ ЮФУ.

Ромм Я.Е, д-р. техн. наук, профессор, заведующий кафедрой информатики ТГПИ.

ISBN 978-5-8327-0268-1   Технологический институт Южного

федерального университета, 2007

         Финаев В.И., Павленко Е.Н.,

Заргарян Е.В., 2007


СОДЕРЖАHИЕ

ВВЕДЕHИЕ……………………………………..……. 6

1. КОНЦЕПЦИЯ МОДЕЛИРОВАНИЯ………...... 9

1.1. Понятие модели…………………………………….. 9

1.2. Концепции определения моделей……………….… 18

 2. КЛАССИЧЕСКИЕ МАТЕМАТИЧЕСКИЕ

МОДЕЛИ…………………………………………..……. 24

2.1. Примеры моделей в виде

дифференциальных уравнений…………………………. 24

2.2. Классические модели в виде

дифференциальных уравнений…………………………. 29

2.3. Инерционные модели…………………………...….. 33

2.4. Модели на основе передаточных функций……..… 38

2.5. Конечные автоматы…………………………………. 40

3. СТОХАСТИЧЕСКИЕ МОДЕЛИ

ОБЪЕКТОВ………………………………………..……. 46

3.1. Математические модели случайных процессов..… 46

3.2. Классификация моделей случайных процессов..… 53

3.3. Модели марковских процессов………………….… 55

 4. ИМИТАЦИЯ СЛУЧАЙНЫХ СОБЫТИЙ….... 59

4.1. Понятие статистического моделирования…...……. 59

4.2. Датчики случайных чисел…………………….…… 60

4.3. Проверочные тесты…………………………..…….. 69

4.4. Имитация случайных событий…………………..… 78

4.5. Имитация непрерывных случайных величин…..… 83

4.6. Имитация марковского процесса……………….…. 93

 5. ОБРАБОТКА РЕЗУЛЬТАТОВ

МОДЕЛИРОВАНИЯ НА ЭВМ………………………. 106

5.1. Выбор числа опытов……………………………….. 106

5.2. Значимость оценки…………………………………. 111

5.3. Формулы и алгоритмы для оценки

результатов моделирования…………………………….. 120

 6. МОДЕЛИРОВАНИЕ ВЕРОЯТНОСТНЫХ

АВТОМАТОВ…………………………………………... 132

6.1. Аналитическое определение

вероятностных автоматов…………………..…………... 132

6.2. Табличное задание функций переходов и выходов 136

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

вероятностных автоматов…………………..…………... 138

7. МОДЕЛИ СИСТЕМ

МАССОВОГО ОБСЛУЖИВАНИЯ……..…………... 147

7.1. Общие сведения…..………………………………..... 147

7.2. Модель входного потока заявок

и времени обслуживания…..…………………….……... 151

7.3. Модель Эрланга…..…………………………….…... 153

7.4. Исследование модели пуассоновского процесса

с помощью производящих функций…..………...……... 156

7.5. Модель для определения времени задержки в виде

интегро-дифференциальных уравнений

Линди-Такача-Севастьянова…..………………………... 158

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

одноканальной СМО…..………………………….……... 162

7.7. Имитационные модели многофазных СМО………. 184

7.8. Имитационные модели многоканальных СМО...… 190

7.9. Алгоритмизация имитационной модели

СМО произвольной структуры…………………………. 211

 8. АЛГОPИТМИЗАЦИЯ ПPОЦЕCCОВ

ФУНКЦИОНИPОВАНИЯ CИCТЕМ………….……. 221

8.1. Моделиpующие алгоpитмы……………………..…. 221

8.2. Пpинципы поcтpоения моделиpующиx алгоpитмов

для cложныx cиcтем………………………………..……. 226

9. УНИФИЦИРОВАННЫЙ

ЯЗЫК МОДЕЛИРОВАНИЯ UML…………..………. 229

9.1. Основные компоненты…………..…………………. 229

9.2. Понятия и компоненты…………..…………………. 231

9.3. Диаграммы вариантов использования ………….… 238

9.4. Диаграммы классов…………………………………. 241

9.5. Типы связей между классами……………………… 249

9.6. Расширения понятия класса в UML………………. 255

9.7. Связи между объектами……………………………. 258

9.8. Диаграммы взаимодействия……………………….. 261

9.9. Диаграммы состояний……………………………… 267

9.10. Диаграммы деятельностей………………………... 278

10. ОБЪЕКТНО-ОРИЕНТИРОВАННОЕ МОДЕЛИРОВАНИЕ………………………..………… 282

10.1. Определение объекта ……………………..……… 282

10.2. Наследование ……………………………………… 295

10.3. Полиморфизм……………………………………… 299

10.4. Типы данных и пакеты………………………….… 303

БИБЛИОГРАФИЧЕСКИЙ СПИСОК…………… 308


ВВЕДЕНИЕ

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

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

Моделирование направлено на решение задач и проведение исследований. При исследовании систем методами системного анализа необходимо построить модель, т.е. реальному объекту ставится в соответствие некоторый математический объект, называемый его моделью. Исследование модели позволяет получить рекомендации относительно поведения реального объекта [1, 2].

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

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

В данном пособии не ставится задача обучения конкретным пакетам прикладных программ для решения задач моделирования, например, Matlab-Simulink-Stateflow, UML, Simulink, Stateflow, Model Vision Studium, Modeliсa, AnyLogiс и других.

Одним из современных направлений компьютерного моделирования является объектно-ориентированное моделирование, которое рассматривают как расширение языка проектирования сложных вычислительных систем — Unified Modeling Language (UML) [4].

В пособии излагается материал по дисциплине «Компьютерное моделирование» согласно государственному образовательному стандарту.

В первом разделе изложены основные концепции моделирования. Рассмотрена классификация моделей, тесно связанная с классификацией систем.

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

В третьем разделе приведено описание стохастических моделей объектов. Рассмотрены способы задания моделей объектов в виде стохастических распределений, классификация стохастических моделей, модели марковских процессов.

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

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

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

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

В восьмом разделе изложены принципы поcтpоения моделиpующиx алгоpитмов процессов функционирования систем.

В девятом и десятом разделах изложены основные сведения об унифицированном языке моделирования UML, об объектно-ориентированном моделировании. За основу взяты материалы работы [4].


1. КОНЦЕПЦИЯ МОДЕЛИРОВАНИЯ

1.1. Понятие модели

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

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

При исследовании систем с применением системного подхода вначале необходимо построить модель, т.е. реальному объекту ставится в соответствие некоторый математический объект, называемый его моделью. Модель (modulus (лат.) - мера) есть объект-заменитель объекта-оригинала [1]. Модель обеспечивает изучение свойств оригинала, а моделирование есть замещение одного объекта другим объектом с целью получения информации о свойствах объекта-оригинала. Теория замещения объектов называется теорией моделирования. Моделирование объекта может быть математическим, в виде построения некоторого макета объекта, натурным и имитационным.

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

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

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

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

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

Имитационное моделирование представляет собой определенную последовательность этапов решения задач:

- изучение реальных систем;

- составление содержательного описания процесса функционирования;

- формулировка цели исследования; выбор основных критериев функционирования;

- разбиение сложной системы на подсистемы;

- построение формализованной схемы процесса функционирования;

- построение математической модели системы;

- планирование эксперимента и сбор исходных данных;

-составление рабочей программы с учетом особенностей машины;

- отладка программы;

- осуществление моделирования;

- обработка результатов;

- выработка рекомендаций.

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

Формализация объекта и разработка адекватных математических моделей – начальная часть работ по моделированию систем самого разного назначения.

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

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

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

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

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

Начало формализации объекта состоит в следующем. Объект функционирует в некоторой среде. Среда воздействует на объект, а объект воздействует на среду. Hа объект могут подаваться управляющие параметры (сигналы). Объект характеризуется вектором состояний Z, называемым еще вектором конструктивных параметров. Формально это отображают схемой, приведенной на рис. 1.1, где Х, Y, F - векторы входных, выходных и возмущающих сигналов.

Рис. 1.1

1.1.2. Классификация моделей. Существует много определений понятия «модель» [1,2,3]. Определим математическую модель как упрощенное отображение существенных сторон реальной системы, выраженное в математической форме и позволяющее математически описать правило (оператор) преобразования входных Х сигналов в выходные Y:

Y=W(Х),       (1.1)

где W – некоторая математическая модель системы.

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

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

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

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

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

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

- аналитически, когда стремятся получить в явном виде зависимости для искомых величин;

- численно, когда нет метода решения уравнения в общем виде, но можно получить результаты при конкретных начальных условиях;

- качественно, когда нет решения в явном виде, но можем найти некоторые свойства решения (оценить устойчивость и т.п.).

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

Модели делятся на стационарные, когда структура и свойства оператора W(t) не изменяются со временем, и на нестационарные.

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

На рис. 1.2 приведен пример изменения выходного сигнала для стационарной и нестационарной системы.

Рис. 1.2

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

Модели делятся на модели безынерционных и инерционных (модели с запаздыванием или с предысторией) систем.

Безынерционные модели соответствуют системам, в которых оператор W определяет зависимость выходных величин от входных в один и тот же момент времени – y(t)=W(Х,t).

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

Y=W(Ztt-1,…,хt-k),      (1.2)

где хt-i=х(t-i) - значение входного сигнала в момент времени t-i. Инерционные модели еще называют моделями с памятью.

Оператор преобразований может содержать параметры, которые обычно неизвестны - Y=W(,Z,Х), где ={1,2,…,k} - вектор параметров.

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

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

 (1.3)

Математическую модель с использованием линейного оператора можно записать в виде Y=WХ. Если условие (1.3) не выполняется, модель называется нелинейной.

Классифицируются динамические модели в соответствии с тем, какие математические операции используются в операторе. Можно выделить: алгебраические, функциональные (типа интеграла свертки), дифференциальные, конечно-разностные модели и др.

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

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

1.2. Концепции определения моделей

Под динамической системой понимается объект, находящийся в каждый момент времени tT в одном из возможных состояний Z и способный переходить во времени из одного состояния в другое под действием внешних и внутренних причин.

При моделировании динамической системы применяют следующие механизмы:

- описание изменения состояний под действием внутренних причин (без вмешательства внешней среды);

- описание приема входного параметра и изменения состояния под действием этого параметра (модель в виде функции перехода);

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

Аргументами входных и выходных параметров системы могут служить время, пространственные координаты, а также некоторые переменные, используемые в преобразованиях Лапласа, Фурье и других.

Входные параметры модели в общем случае могут быть заданы в виде вектора ={х12,…,хm}, хiХi,, (), где Хi - заданные дискретные или непрерывные множества. Прямое произведение вида Х=Х1Х2Хm называется пространством входных параметров, а вектор  представляет собой точку пространства Х.

Отображение Х=L(t), сопоставляющее каждому моменту времени t некоторый параметр Х, называется входным процессом L(t).

Вектор выходных параметров  множеству выходных параметров. Выходной параметр, выдаваемый системой в момент времени tT, обозначим . Если выходной сигнал  описывается набором характеристик y1,y2,…,yr, таких, что yjYj, (), где Yj - заданные множества, то прямое произведение Y=Y1Y2Yr называется пространством выходных параметров. По аналогии с входным процессом определяется понятие выходного процесса Y=M(t).

В теории управления выходные параметры называются фазовыми координатами (переменными состояния).

Состояние системы определяется как совокупность состояний элементов. Состояние системы  описывается некоторым набором характеристик zkZk, (), где Zk - заданные множества, а пространство состояний Z определяется как прямое произведение Z=Z1Z2Zn.

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

Рис. 1.3

Выходные сигналы системы и дополнительные воздействия, которым соответствует r-мерный вектор дополнительных сигналов (связанных также с ошибками измерения) ={1,2,…,r}, являются входными сигналами для измерительного устройства. Наблюдаемый вектор состояний измерительной системы (вектор откликов) записывается в виде V={v1,v2,…,vr}. Математическая модель измерительного устройства имеет вид

V=B(Y),       (1.4)

где B(Y) - некоторый оператор, преобразующий сигналы Y и на входе измерительного устройства в сигналы-отклики V.

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

Соответствие [4]  это способ (закон) сопоставления элементов хХ с элементами yY так, что имеется возможность образования пар (двоек) (х,y), причем для каждого элемента хХ возможно указать элемент yY, с которым сопоставляется элемент х. В сопоставлении могут участвовать не все элементы Х и Y. Для задания соответствия необходимо указать:

- множество Х, элементы которого сопоставляются с элементами другого множества;

- множество Y, с элементами которого сопоставляются элементы множества Х;

- множество QХY, определяющее закон, согласно которому осуществляется соответствие, т.е. перечисляющее все пары (х,y), участвующие в сопоставлении.

Соответствие, обозначаемое через q, представляет собой тройку множеств q=(Х,Y,Q), где Х – область отправления соответствия, – область прибытия соответствия, Q – график соответствия, QХY. Очевидно, что проекция Пр1QХ, а Пр2QY, причем множество Пр1Q называется областью определения соответствия, а проекция Пр2Q – областью значений соответствия. Способы задания соответствий следующие.

 При теоретико-множественном задании определяют множества Х={х12,…,хn}, Y={y1,y2,…,ym} и график Q={(хi,yj)}, хХ, yY , .

 При матричном способе задания соответствие задается в виде матрицы инцидентности RQ, которая имеет вид прямоугольной таблицы размером nm. Элементы хiХ соответствуют строкам матрицы RQ, а элементы yjY соответствуют столбцам. На пересечении хi строки и yj столбца ставится элемент rij=1, если элемент i,yj)Q, и rij=0, если i,yj)Q.

 При графическом способе соответствие задается в виде рисунка (см. рис. 1.4.), на котором элементы хiХ – кружки одной линии, элементы yjY – кружки другой линии, а каждая двойка i,yj)Q обозначается стрелкой, идущей от кружка хi к кружку yj. Такое представление называется графиком.

Если сопоставлять элементы yY элементам множества Х, то получим соответствие q-1=(Y,Х,Q-1), обратное соответствию q (инверсия соответствия q).

Х={х1234}, Y={y1,y2,y3}, Q={(х1,y1), (х1,y2), (х2,y1), (х2,y2), (х3,y2), (х4,y3)}.

Рис. 1.4

Исходя из приведенных выше определений множеств входных параметров Х=Х1Х2Хm, выходных параметров Y=Y1Y2Yr, состояний Z=Z1Z2Zn определим задание моделей функций переходов и выходов, как соответствий. Модель системы в виде функции переходов задана соответствием

fП=(Х1Х2Хm,Z1Z2Zn,FП).    (1.5)

Данная модель устанавливает соответствие fП между каждым элементом ={х12,…,хm}Х1Х2Хm и элементом ={z1,z2,…,zn}Z1Z2Zn. FП – график соответствия fП. Модель системы в виде функции переходов может быть задана также в виде

fП={(Х1Х2Хm),(Z1Z2Zn),

(Z1Z2Zn),FП},      (1.6)

т.е. модель устанавливает соответствие fП между каждым элементом (,)[(Х1Х2Хm)(Z1Z2Zn)] и элементом ={z1,z2,…,zn}Z1Z2Zn. Модель системы в виде функции выходов задана соответствием

fВ={[(Х1Х2Хm),(Z1Z2Zn)],

(Y1Y2Yr),FВ}.     (1.7)

Модель устанавливает соответствие fВ между каждым элементом (,) из множества [(Х1Х2Хm), (Z1Z2Zn)] и элементом ={y1,y2,…,yr}Y1Y2 Yr. FВ – график соответствия fВ.

Модель системы в виде функции выходов может быть задана еще и в таком виде:

fВ={[(Х1Х2Хm)(Z1Z2Zn)],[Z1Z2Zn],

(Y1Y2Yr),FВ),      (1.8)

т.е. модель в данном случае устанавливает соответствие fВ между каждым элементом {(,),} из множества {[(Х1Х2Хm), (Z1Z2Zn)], [Z1Z2Zn]} и элементом ={y1,y2,…,yr}Y1Y2 Yr.

Модель системы в виде функции переходов может быть записана еще в следующем виде:

  (1.9)

или в виде

. (1.10)

Модель системы в виде функции выходов может быть задана и в таком виде:

 (1.11)

или в виде

.  (1.12)


2. КЛАССИЧЕСКИЕ МАТЕМАТИЧЕСКИЕ МОДЕЛИ

2.1. Примеры моделей в виде дифференциальных уравнений

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

Классическая теория автоматического управления применяет модели объектов в виде дифференциальных уравнений и позволяет решать задачи управления, исследуя дифференциальные уравнения [5]. Вид уравнений может быть очень сложным и для решения и исследования моделей, заданных дифференциальными уравнениями, применяется компьютерное моделирование. Применяют специальные пакеты прикладных программ, например MatLab, Omola, Dymola, Dymosim, Model Vision Studium, язык UML и другие.

Колебательный контур, известный в электротехнике, представляют в виде схемы, приведенной на рис. 2.1. Определены параметры колебательного контура:

С – емкость конденсатора;

L – индуктивность катушки;

UС(t) – напряжение на конденсаторе;

IL(t) – ток в катушке;

U(t)ИСТ – напряжение внешнего источника.

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

Рис. 2.1

В соответствии с законом Кирхгофа запишем следующие соотношения:

,   .

Введем координаты z1=UС, обозначим UИСТ/L=х(t) и получим систему дифференциальных уравнений:

,   .   (2.1)

Если UИСТ=0, то х(t)=0 и система (2.1) описывает свободные колебания. Рассматривая х(t) как сигнал управления, получим описание динамики колебаний в каждый момент времени t. Решая систему (2.1), можно получить аналитический вид функций z1(t) и z2(t).

Если рассматривать U(t)ИСТ как единичное ступенчатое воздействие, R - внутреннее сопротивление источника, то можно показать, что напряжение на конденсаторе С будет изменяться по формуле [6]:

,    (2.2)

где , , .

Дальнейшее исследование модели связано с разработкой программы, которая будет рассчитывать и выводить график функции UС(t) в зависимости от параметров R, L и С, тем самым позволяя исследовать время переходного процесса (время затухания колебаний) в данном контуре.

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

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

Переходим к предельному выражению

и получаем модель роста популяций микроорганизмов в виде дифференциального уравнения (общий вид):

.       (2.3)

Решая дифференциальное уравнение (2.3), можно получить аналитическое уравнение роста популяций и провести исследования. При начальных условиях t=0, E(t=0)=E0 получим решение модели роста популяций в виде аналитического выражения

E(t)=E0ekt.         (2.4)

Вид уравнения (2.4) показан на рис. 2.2.

Рис. 2.2

Если при t=0 начальное число популяций микроорганизмов E=E0, то можно определить время Т, за которое число особей удвоится по формуле

2E0=E0ekt, 2=ekT, T=(1/k)ln2.

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

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

Пусть m1 - число боевых единиц красных; m2 - число боевых единиц синих, сохранившихся непораженными к моменту времени t; λ1 - средняя скорострельность для одной боевой единицы красных; λ2 - средняя скорострельность для одной боевой единицы синих. Цели поражаются с вероятностью р1 - красными и вероятностью р2 - синими. Рассмотрим модель, отображающую динамику боя.

Интенсивности успешных выстрелов определятся как

L11р1, L22р2.

Число выведенных боевых единиц красных m1 за время t составит λ2р2tm2, а число выведенных из строя боевых единиц синих m2 за время t составит λ1р1tm1, Тогда

m1=-λ2р2tm2,  m2=-λ1р1tm1.    (2.4)

Уравнения (2.4) - модель динамики боя в частных приращениях. От уравнения (2.4) осуществим переход к дифференциальным уравнениям.

Разделив правую и левую части на t, получим

,    .

Взяв пределы при t0, получим дифференциальные уравнения, моделирующие динамику боя

,       (2.5)

Уравнения (2.5) называются уравнениями Ланчестера.

2.1.4. Модель движения ракеты. Движение ракеты, запускаемой в космос, описывается её координатами х и y, проекциями вектора скорости V на координатные оси VХ и VY. Пусть m - масса ракеты; u величина тяги; - угол между направлением тяги и осью ; f(u) - секундный расход массы. Рассмотрим построение модели, отображающей динамику полета.

Проекции скоростей являются производными от движения по координатам, следовательно

,  .

В соответствии с уравнением Ньютона запишем:

  .

Расход массы определится уравнением

.

Таким образом, моделью движения ракеты является система уравнений:

при начальных условиях х(t0)=х0, y(t0)=y0, m(t0)=m0, Vх(t0)=Vх0, Vy(t0)=Vy0.

Управление траекторией ракеты осуществляется за счет регулирования величины и направления силы тяги двигателя, U и - управляющие параметры.

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

2.2. Классические модели в виде дифференциальных уравнений

Дифференциальные уравнения описывают процесс перехода динамической системы из одного состояния в другое и изменение выходного параметра. Могут рассматриваться системы, в которых моделируют только изменение состояний или только изменение выходного параметра. Могут рассматриваться системы, в которых моделируют изменение и состояний, и выходного параметра. Математические схемы такого вида отражают динамику изучаемой системы и носят название D-схем от английского слова dynamic (динамика).

Пусть входные параметры (сигналы, координаты и прочее) заданы множеством Х(t)={х1(t),х2(t),...,хm(t)}, а выходные параметры заданы множеством Y(t)={y1(t),y2(t),..., yr(t)}. Модель динамической системы, определяемая обыкновенными дифференциальными уравнениями в общем виде, задается следующим образом.

Задают дифференциальные уравнения, определяющие движение системы в пространстве состояний

.  (2.6)

Каждое i-е дифференциальное уравнение задается в общем виде функцией fi, зависящей от времени t, компонент вектора состояний Z={z1(t),z2(t),…,zn(t)} и компонент вектора входных параметров Х(t)={х1(t),х2(t),...,хm(t)}. Задают соотношения, определяющие изменение выходных параметров

.  (2.7)

Для решения дифференциальных уравнений системы (2.6), определения изменения во времени выходных параметров необходимо для момента t(0)=t0 задать начальные состояния , 

а также функции, определяющие изменения во времени компонент вектора входных параметров Х(t) на полуинтервале (t0,t]:

.

Если для каждого уравнения системы (2.6) выполнены условия существования и единственности решений, то эти решения в общем случае имеют вид

.  (2.8)

Обозначим решения системы дифференциальных уравнений (2.6), проходящие в момент времени t0 через точку , символом F. Тогда модель в виде функции переходов для динамической системы будет задана в общем виде уравнением

.      (2.9)

Эта функция каждому набору  ставит в соответствие то состояние Z(t), в которое переходит система за время перехода t-t0 из фазы (t0,Z0) под действием входных параметров .

Модель динамической системы в виде функции выходов в общем виде будет определена уравнением

   (2.10)

в котором оператор G каждому набору  сопоставляет выходной сигнал yt=y(t).

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

Если модель предназначена для описания изменения состояния z(t) динамической системы, то модель в виде обыкновенного линейного дифференциального уравнения q-го порядка с постоянными коэффициентами и правой частью, выраженной через производные от управляющих функций, задается в следующем виде:

  (2.11)

Если применить оператор дифференцирования , то с учетом аддитивной ошибки v(t) уравнение (2.11) запишется в виде

z(р)=-1(р)(р)х(р)+v(р),

где -1(р)=рq-1рq-1-2рq-2-…-q,  (р)=0рr+1рr-1 + … + r.

Модели в виде многомерных дифференциальных уравнений в форме Коши находят наибольшее применение. Они описываются системами обыкновенных дифференциальных уравнений первого порядка в форме Коши, т.е. разрешенными относительно первых производных. Для стационарной линейной системы, параметры которой изменяются непрерывно во времени, модель в общей форме имеет следующий вид:

.   (2.12)

В уравнении (2.12):  Z={z1(t),z2(t),…,zn(t)} - вектор состояний; Х(t)={х1(t),х2(t),...,хm(t)} – вектор входных параметров; Y(t)={y1(t),y2(t),..., yr(t)} – вектор выходных параметров; W={w1(t),w2(t),…,wn(t)} - вектор шума системы;  - транспонированный вектор производных от переменных состояния; матрицы Ф, G, H и Г имеют размерности, зависящие от размерностей векторов Z, Х(t), Y(t), W. Коэффициенты матриц Ф, G, H и Г имеют смысл коэффициентов передачи, для стационарной системы не зависят от времени и подлежат оцениванию. Параметры могут входить и в начальное условие, которое необходимо добавить для решения первого уравнения (2.12).

Модель для нестационарной линейной непрерывной системы отличается от (2.12) тем, что матрицы Ф, G, H и Г будут зависеть от времени.

Непрерывная нелинейная система может быть описана моделью

(2.13)

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

Компьютерное моделирование систем, описываемых многомерными дифференциальными уравнениями в форме Коши, осуществляется с применением пакетов программ. Широко используется подсистема Simulink пакета MatLab. При моделировании определяется вид дифференциального уравнения, задаются начальные условия. Результаты решения отображаются визуально в виде цифровых данных, а также в виде графических данных.

2.3. Инерционные модели

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

2.3.1. Дифференциальные уравнения с запаздывающими аргументами. В общем случае дифференциальные уравнения n-го порядка с запаздывающим аргументом имеют вид

.     (2.14)

Дифференциальное уравнение (2.14) может быть сведено к системе дифференциальных уравнений первого порядка. Введем обозначения: z=z1; ; и т.д. Тогда дифференциальное уравнение (2.14) запишем в следующем виде:

.

Из рассмотрения даже простейшего дифференциального уравнения вида

,    (2.15)

где >0, =сonst, сложно понять, какие начальные условия надо задать, чтобы определить решение z(t) для t>t0.

Решение дифференциального уравнения (2.15) определяется из интегрального уравнения

 (2.16)

Решение уравнения (2.16) осуществляется по следующему алгоритму.

Следует задать начальное значение для точки t0 z0=z(t0) и функцию z(t) в полуинтервале t0-t<t0 ([t0-, t0)). Функцию z(t)=W(t) называют начальной функцией t[t0-,t0). При таких условиях можно получить либо аналитическое решение уравнения (2.16), либо получить решение для любого >t0 с применением методов вычислительной математике и компьютерного моделирования. Алгоритм решения уравнения (2.16) представляет собой следующую последовательность действий.

После задания начальных условий следует определить непрерывное решение z(t) для t>t0, при условии, что z(t)=W(t) для t[t0-,t0). Если функции f и W непрерывны и первая из них удовлетворяет условию Липшица по z, то искомое решение существует и единственно.

Зная W(t) для t0-t<t0, найдем z(t) для t0t<t0+. Примем это z(t) за начальную функцию W(t) для t0t<t0+. Определим z(t) для t0+t<t0+2 и т.д.

При поиске решения применен метод последовательного интегрирования, сущность которого показана на рис. 2.3.

Дифференциальные уравнения с запаздывающим аргументом применяются для составления моделей динамической системы с последствием, т.е. систем, для определения состояний z(t) которых при t>t0 недостаточно задать z0=z(t0).

Рис. 2.3

2.3.2. Модели в виде сумм и интегралов свертки. Если динамическая система функционирует в дискретные моменты времени, то ее модель может быть описана в виде суммы свертки. Математические модели, выражаемые суммой свертки или интеграла свертки, задаются следующим образом. Для однооткликовой стационарной динамической системы, на вход которой действует управляющая функция х(t), а наблюдения над входом и выходом производятся только в дискретные моменты времени с интервалом квантования t, математическая модель может быть выражена с помощью суммы свертки:

Зададим масштаб t=1, получим

    (2.17)

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

Если линейная динамическая система нестационарная, то в выражения (2.17) нельзя применять импульсную характеристику системы вида h(k-i). Для этого случая математическая модель примет вид

где h(k,i) - реакция системы в момент k на единичный импульс в момент i.

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

- для линейной системы

; (2.18)

- для нестационарной системы

Модель представлена в виде функционала с аддитивной ошибкой v(t). Интеграл называется интегралом свертки, или интегралом Дюамеля. Для определения импульсной характеристики (весовой функции) используется (для стационарных систем) представление весовой функции в форме Релея-Ритца путем разложения функций в ряд по системе известных ортогональных функций

    (2.19)

где Фi(t) - функции системы ортогональных функций при значениях t, принадлежащих отрезку ортогональности [t1,t2]. Это позволяет сделать модель параметрической, которая содержит ограниченное число параметров i, подлежащих определению. Коэффициенты i называют еще спектром разложения в ряд базисных функций.

К системе базисных функций предъявляются следующие требования: для любой весовой функции ряд (2.19) должен сходиться; Фi(t) должна иметь простую аналитическую форму; i должны вычисляться аналитически просто.

Условие ортогональности базисных функций имеет вид

,    (2.20)

где число сi называют нормой базисной функции Фi(t).

Каждую базисную функцию можно нормировать по ее норме, причем нормированная функция имеет вид

.

Система (2.20) примет вид

,  (2.21)

где ij - символ Кронекера.

Для определения i умножим правую и левую части уравнения (2.19) на Фi(t) и проинтегрируем обе части на отрезке ортогональности

.

При k=i правый интеграл равен единице, тогда

.      (2.22)

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

2.4. Модели на основе передаточных функций

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

Дискретное Z-преобразование решетчатой функции v(n)=vn задается в следующем виде

.

Применяя одностороннее Z-преобразование к левой и правой части этого выражения, получаем

/     (2.23)

Z-преобразование однозначно связано с дискретным преобразованием Лапласа. Взаимосвязь комплексной переменной z и комплексной переменной преобразования Лапласа выражается соотношением z=es, которое используется для перехода от дискретного  преобразования  Лапласа к Z-преобразованию и наоборот.

Модель импульсной системы (2.23) устанавливает связь между Z-преобразованием  отклика z(k) выходного сигнала и Z-преобразованием  входного сигнала х(k).  - передаточная функция импульсной системы (дискретная передаточная функция), являющаяся Z-преобразованием импульсной характеристики h(k).  - Z-преобразование случайной составляющей v(k).

Если применять преобразование Лапласа к обеим частям модели (2.18) для непрерывной однооткликовой системы, то можно записать z(s)=h(s)х(s)+v(s). В этом уравнении z(s), h(s), х(s), v(s) - преобразования Лапласа соответственно от z(t), h(t), х(t), v(t); h(s) - передаточная функция непрерывной системы, представляющая собой преобразование Лапласа от импульсной характеристики. Преобразование Лапласа имеет вид

.

Применяя к обеим частям уравнения (2.18) дискретное преобразование Фурье, получим z(jw)=h(jw)х(jw)+v(jw), где z(jw), х(jw), v(jw) - преобразования Фурье соответственно от отклика, входного сигнала и помехи, h(jw) - частная характеристика системы, которая есть не что иное, как преобразование Фурье от импульсной характеристики. Преобразование Фурье имеет вид

.

В рассмотренных моделях, использующих преобразования по Лапласу и Фурье, в роли аргументов выступают соответствующие параметры преобразований z, s, j. Все модели линейны по входным сигналам, но, как правило, нелинейные по параметрам.

2.5. Конечные автоматы

Для моделирования динамических систем, функционирующих в дискретном времени, применяется аппарат конечных автоматов [7]. Теория конечных автоматов и их модели используются при синтезе и анализе вычислительных устройств, дискретных устройств автоматики. Конечный автомат функционирует в дискретные моменты времени t, причем в каждый момент ti автомат находится в одном из возможных состояний z(ti), принадлежащем множеству состояний автомата Z. Математические модели в виде конечного автомата получили название F-схем от английского finite automata – конечный автомат.

В каждый момент ti (i=1,2,...) на вход конечного автомата поступает входной параметр - одна из букв х(ti) входного алфавита Х, а на выходе существует выходной параметр y(t) - буква выходного алфавита Y.

Автомат формально определен набором

A=<Х,Z,Y,z0,,>,

где Х={х12,...,хm} - множество входных параметров; Z={z1,z2,...,zn} - множество состояний; Y={y1,y2,...,yr} - множество выходных параметров. Элементы множества Х, Z, Y называют входным, внутренним и выходным алфавитом. При поступлении параметра х состояние конечного автомата изменяется в соответствии с одношаговой функцией переходов, например:

z(t)=[z(t-1),х(t)] или z(t)=[z(t),х(t)],

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

y(t)=[z(t),х(t)]; y(t)=[z(t-1),z(t) х(t)]; y(t)=[z(t-1), х(t)].

Функции переходов и выходов могут быть заданы теоретико-множественным способом, табличным способом и в виде графов. Рассмотрим пример задания конечного автомата.

Пусть Х={х123}, Z={z1,z2,z3,z4}, Y={y1,y2,y3,y4}.

Функция переходов для данного автомата задана в виде табл. 2.1, а функция выходов вида y(t)=[z(t),х(t)] – в виде табл. 2.2.

                Таблица 2.1

Функция переходов

Х

Z

z1

Z2

z3

z4

х1

z1

Z2

z3

z4

х2

z2

Z3

z4

z1

х3

z4

Z1

z2

z3

                Таблица 2.2

Функция выходов

Х

Z

z1

z2

z3

z4

х1

y1

y4

y3

y2

х2

y2

y1

y4

y3

х3

y3

y2

y1

y4

При задании функции переходов на пересечении i–й строки и j–го столбца указывается состояние zk, в которое переходит автомат при подаче входного параметра хi в такте времени t, при условии, что в такте времени t-1 автомат находился в состоянии zj.

На рис. 2.4 приведено графическое задание функции переходов рассматриваемого автомата.

Рис. 2.4

Автомат в процессе своей работы реализует отображение множества слов (последовательность параметров) входного алфавита Х на множество слов выходного алфавита Y. Если на вход конечного автомата, установленного в начальное состояние z0, подать последовательность букв входного алфавита х(t0), х(t1), х(t2),…, то на выходе автомата будут последовательно появляться буквы выходного алфавита y(t0), y(t1), y(t2),…

В зависимости от способа заданий функций переходов и выходов, автоматы подразделяются на автоматы первого и второго рода. Для автомата первого рода, называемого автоматом Мили, функция переходов имеет вид z(t)=[z(t),х(t)], а функция выходов - y(t)=[z(t),х(t)].

Для автомата второго рода функция переходов имеет вид: z(t)=[z(t),х(t)], а функция выходов - y(t)=[z(t-1),х(t)].

Автомат второго рода, функция выходов которого определяется его состоянием - y(t)=[z(t)], называется автоматом Мура.

По числу состояний различают конечные автоматы с памятью и без памяти. Автоматы с памятью имеют более одного состояния, а автоматы без памяти обладают лишь одним состоянием. Автоматы без памяти называют комбинационными или логическими схемами. Функция выходов такого автомата - y(t)=[х(t)].

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

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

Пример 1. Автомат для продажи билетов в автобусах принимает монеты достоинством в 1,2 и 5 рублей и выдает билеты стоимостью 5 рублей. Рассмотрим конечный автомат Мили с множеством состояний Z=(0,1,2,3), входным алфавитом Х=(1,2,5) и выходным алфавитом Y=(0,1), где 0 соответствует ситуации «билет не выдается», а 1 – ситуации «билет выдается».

Функция переходов (t) определяется соотношением

Z(t)=(z(t-1)+х(t))mod5,

а функция выходов (t) – соотношением

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

Изделия i-й номенклатуры хранятся в i-м стеллаже (i=1,2,…,n). Содержание стеллажей изменяется в моменты времени поступления на склад новых изделий потребителям. Такое хранилище можно представить в виде конечного автомата Мура. В качестве состояний выберем n-мерный вектор Z=(Z1,Z2,…,Zn), где Zi - число изделий на i-м стеллаже.

Выходной сигнал – (n+1) – мерный вектор Х=(Х12,…, Хn,), где Хi - число изделий i-й номенклатуры, поступивший на склад или выданный потребителю. При поступлении изделий на склад =1, а при выдаче изделий потребителю = -1. Выходной сигнал представляет собой n–мерный вектор Y=(Y1,Y2,…,Yn), для которого Yi(t)=Zi(t) – информация о состоянии стеллажей.

Функция переходов определена соотношением

Zi(t)=Zi(t-1)+х(t),

а функция выходов определена соотношением

Yi(t)=Zi(t).

На рис. 2.5 приведен алгоритм программы моделирования конечного автомата, функция выходов которого имеет задание z(t)=[z(t-1),х(t)], а функция переходов - y(t)=[z(t),х(t)]. Функция переходов задается в подпрограмме WWOD в виде массива FР(i,j), значения которого определены индексом k состояния zk, в которое переходит автомат при подаче входного параметра хi в такте времени t, при условии, что в такте времени t-1 автомат находился в состоянии zj.

Рис. 2.5

Функция выходов задается в подпрограмме WWOD в виде массива FW(i,k), значения которого определены индексом w выходного параметра (буквы) yw, который будет на выходе автомат при подаче входного параметра хi в такте времени t, при условии, что в этом такте времени t автомат находится в состоянии zk. В подпрограмме WWOD происходит инициализация автомата, т.е. определение z0=zk.  В блоке 3 алгоритма происходит изменение индекса состояния с учетом задержки на один такт.


3. СТОХАСТИЧЕСКИЕ МОДЕЛИ ОБЪЕКТОВ

3.1. Математические модели случайных процессов

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

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

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

3.1.1. Распределение случайных событий. Массовые явления или процессы характеризуются многократным повторением при постоянных условиях некоторых опытов (операций и прочее). Абстрагируясь от специальных свойств этих опытов, в теории вероятностей вводится понятие испытания (опыта). Испытанием называется осуществление определенного комплекса условий, который может быть воспроизведен сколь угодно большое число раз. Явления, происходящие при реализации этого комплекса условий (в результате испытания), называются событиями [8].

Положительное число в отрезке [0,1], представляющее собой количественную меру возможности реализации случайного события в испытании, называется его вероятностью. Вероятность появления события А обозначают символом Р(А), причем 0Р(А)1. Вероятность понимается как идеальная мера возможности появления события.

Случайная величина рассматривается как функция, аргументом которой служит элементарное случайное событие. Дискретной случайной величиной называется такая, которая может принимать конечное или бесконечное счетное множество значений, например возможны значения x1, x2, …, xn, … Для каждого события xi определены вероятности P(xi). Распределение вероятностей дискретной случайной величины, представленное на рис. 3.1, рассматривают как точечное распределение вероятностей.

Рис. 3.1

При непрерывном распределении случайной величины вероятности распределены сплошной полосой по всей оси x или по некоторым ее участкам с определенной плотностью.

Распределение вероятностей носит название теоретического распределения случайной величины.

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

.     (3.1)

Пример задания интегральной функции распределения вероятностей приведен на рис. 3.2.

Рис. 3.2

Дифференциальная функция распределения вероятностей (плотность распределения вероятностей) определяет вероятность того, что случайная величина X меньше значения x

.  (3.2)

Пример задания дифференциальной функции распределения вероятностей приведен на рис. 3.3.

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

Рис. 3.3

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

Определение. Моделью случайного процесса называют случайную функцию X(), заданную на множестве , принимающую действительные значения и описываемую семейством распределений [9]:

, i, i=1,2,...,n, n=1,2,...,

которое удовлетворяет условиям согласованности

,

=,

где i1, i2,…, in, - любая перестановка индексов 1, 2,..., n.

Набор функций  называется конечномерными распределениями случайной функции или интегральной функции распределения вероятностей многомерной случайной величины. При n=1 получим одномерное распределение (3.1). Модель многомерного распределения необходима для моделирования многопараметрической случайной величины.

При решении многих задач моделирования приходится оперировать с несколькими случайными функциями. Для того чтобы над ними производить математические операции, недостаточно, чтобы каждая из этих случайных функций была задана в отдельности. Последовательность функций X1(), X2(),…, Xn() возможно заменить векторной функцией (), компонентами которой служат случайные функции Xi(), (i=1,2,…,n).

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

Если  - плотность функций распределения , то

=

=.

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

.

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

X1(),X2(), …, Xn(), i0 >, i=1,n, n=1,2,...,

которая определяется формулой

,

где M - символ математического ожидания, u1,u2,...,uk - вещественные числа.

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

.

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

Момент k–го порядка дискретной случайной величины определяется по формуле

.

Для непрерывной случайной величины моментная функция k–го порядка определяется по формуле

.

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

Определение. Модель случайной функции X(i), i в виде моментной функции задается отношением

,

если математическое ожидание в правой части равенства имеет смысл при всех i, i=1,n. Величина q=j1+j2+...+jn называется порядком моментной функции.

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

при u1=u1=…=un=0.

Кроме моментных функций в качестве моделей часто рассматривают центральные моменты функции. Центрированной случайной величиной называется случайная величина . Для непрерывной случайной величины центральная моментная функция k–го порядка определяется по формуле

.

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

,

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

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

m()=m1(1)=MX(),

R1(1,2)=m1(1,2)=M{[X(1)–m(2)][X(2)–m(2)]}.

Функции m() называются средним значением или математическим ожиданием, а R1(1,2) - корреляционной функцией. При 1=2= корреляционная функция дает дисперсию () величины (), R1(1,2)=2().

Величину

называют коэффициентом корреляции случайных величин X(1) и X(2).

3.2. Классификация моделей случайных процессов

Случайные процессы делятся на следующие широкие классы: гауссовы процессы; процессы с независимыми приращениями; стационарные в широком смысле; марковские процессы.

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

Определение. Случайный вектор X=(X1,X,...,Xn) имеет гауссово (нормальное) распределение, если характеристическая функция распределения представима в виде

(u)=M{eхр[j(u,x)]}=eхр[j(m,u)-0,5R(u,u)],

где m=(m1,m2,…,mn), u=(u1,u2,…,un) - векторы, R - неотрицательно-определенная вещественная симметричная матрица, R=||rij||, i,j=1,n. Здесь (,) обозначает скалярное произведение векторов  и , так, что

,  .

3.2.2. Модель процессов с независимыми приращениями. Пусть T - конечный отрезок T=[0,a] или T=[0,].

Определение. Случайный процесс {X(t), tT} со значениями в евклидовом пространстве Rn называется процессом с независимыми приращениями, если для любых n, таких, что 0<t1<t2<...<tn, случайные векторы X(0), X(t1)-X(0),...,X(tn)-X(tn-1) - взаимно независимы.

Вектор X(0) называется начальным состоянием (значением) процесса, а его распределение - начальным распределением процесса. Чтобы задать процесс с независимыми приращениями в широком смысле, достаточно задать начальное распределение Р0(B)=Р{X(0)B} и набор распределений Р(t,h,B) - распределений вектора Р{X(t+h)-X(t)}B.

Процесс с независимыми приращениями называется однородным, если распределения вектора X(t+h)-X(t) не зависят от t, Р(t,h,B)=Р(h,B).

3.2.3. Модель процессов, стационарных в широком смысле. Стационарные процессы - это такие процессы, теоретико-вероятностные характеристики которых не изменяются со временем. Пусть T=[0,a] или T=[0,).

Определение. Модель случайного процесса (в широком смысле) {X(t), tT} со значениями в Rn называется стационарной, если для любого n и любых t1,t2,...,tт, таких, что tk+tT, (k=1,n), совместное распределение случайных векторов, описывающих случайный процесс X(t1+t),...,X(tn+t), не зависит от t.

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

Определение. Случайный процесс X(t), t>0 со значениями в пространстве Rn называют процессом, стационарным в широком смысле, если M[X(t)]2< и M[X(t)]=monst, M[X(t)-m][X(s)-m]=R(t-s), (t>s), где R(t) - непрерывная матричная функция.

Функцию R(t) называют корреляционной (матричной) функцией процесса X(t). В качестве примера стационарных в широком смысле процессов можно рассмотреть колебания со случайными параметрами.

3.3. Модели марковских процессов

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

Представим себе систему, которая может находиться в разных состояниях. Возможные состояния образуют множество Х, называемое фазовым пространством. Пусть система эволюционирует во времени. Состояние системы в момент времени t обозначим через хt. Если хtB, где BХ, то будем говорить, что система в момент времени t находится во множестве B. Предположим, что эволюция системы носит стохастический характер, т.е. состояние системы в момент времени t не определяется однозначно через состояние системы в моменты времени s, предшествующие t, где s<t, а является случайным и описывается теоретико-вероятностными законами.

Пусть Р(s,х,t,B) - вероятность события хtB (s<t), при условии, что хs. Функцию Р(s,х,t,B) называют вероятностью перехода рассматриваемой системы. Под системой без последействия понимают систему, для которой вероятность попадания в момент времени t во множество B, при полностью известном движении системы до момента времени s (s<t), по-прежнему равна Р(s,х,t,B) и, таким образом, зависит только от состояния системы в последний момент времени.

Обозначим через Р(s,х,u,y,t,B) условную вероятность события хtB при гипотезах хs=х, хu=y (s<u<t). В соответствии с общими свойствами условных вероятностей имеет место равенство

. (3.3)

Для системы без последствия естественно предположить, что

Р(s,х,u,y,t,B)=Р(u,y,t,B).

Тогда равенство (3.3) примет вид

.  (3.4)

Соотношение (3.4) называется уравнением КолмогороваЧепмена. Это уравнение определяет модель марковского процесса.

Пусть {Х,B}-некоторое измеримое пространство. Функцию Р(х,B), хХ, BB, удовлетворяющую условиям:

а) Р(х,B) при фиксированном х является мерой на B и Р(х,Х)=1;

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

Пусть I - некоторый конечный или бесконечный полуинтервал (отрезок). Семейство стохастических ядер st(х,B)=Р(s,х,t,B), s<t, (s,t)II}, удовлетворяющих уравнению Колмогорова-Чепмена (3.4), будем называть марковским семейством стохастических ядер.

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

- измеримое пространство {х, B};

- полуинтервал I (отрезок) действительной оси;

- марковское семейство стохастических ядер st(х,B), s<t, (s,t)II}.

Семейство ядер Рst(х,B)=Р(s,х,t,B) называют вероятностью перехода марковского процесса, пространство {х,B} - фазовым пространством системы, точка множества I интерпретируется как моменты времени, а величина Рst(х,B)=Р(s,х,t,B) - как условная вероятность того, что система в момент времени t окажется во множестве B, если в момент времени s она находилась в точке х фазового пространства (s<t).

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

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

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

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

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


4. ИМИТАЦИЯ СЛУЧАЙНЫХ СОБЫТИЙ

4.1. Понятие статистического моделирования

При определении методов статистического моделирования применяют название «метод Монте-Карло». Определение, которое характеризует этот метод достаточно точно и полно, не существует. Известно, что этот метод всегда связан со случайными испытаниями. Все расчеты по методу Монте-Карло заключаются в случайной выборке из некоторой «генеральной совокупности» в соответствии с определенными вероятностными законами. Метод Монте-Карло появился в 1949 году. В этом году вышла в свет статья «The Monte Carlo method» американских авторов Metropolis N. и Ulam S. Создателями этого метода считают математиков Дж. Неймана и С. Улама [10].

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

Название «Монте-Карло» происходит от города Монте-Карло в княжестве Монако, знаменитого своим игорным домом, т.к. рулетка служила простейшим механическим прибором для получения случайных величин.

Первоначально метод успешно применялся для решения задач ядерной физики - переноса нейтронов, расчета критичности ядерных реакторов, задач защиты от излучений и других. Широкое применение данного метода при решении различных задач началось с развития вычислительной техники, особенно при появлении персональных ЭВМ.

Метод Монте-Карло имеет две отличительные особенности. Первая заключается в простоте структуры вычислительного алгоритма - многократного повторения одного случайного испытания. Так как испытания являются независимыми опытами, то результаты всех испытаний усредняются. В связи с этим метод иногда называют методом статистических испытаний. Вторая особенность состоит в зависимости точности метода от числа испытаний N. Метод достаточно эффективен при решении задач, не требующих высокой точности, например, в пределах 5-10%. Метод статистического моделирования позволяет решать как детерминированные, так и вероятностные задачи [11].

4.2. Датчики случайных чисел

Для имитации случайных событий необходим некоторый эталон, т.е. то, с чем можно что-то сравнить. Известно, что наука существует там, где есть измерения. Отсутствие измерений приводит к схоластике, лженаукам. В мире существуют эталоны для измерения времени, расстояния, веса и прочее. Такой же эталон должен существовать и в ЭВМ для «измерения» вероятности случайного события. Начало созданию подобного эталона было положено в шестидесятых – семидесятых годах прошлого столетия. Создавались специальные устройства – датчики (генераторы) случайных чисел или случайных величин, авторы которых получали авторские свидетельства на эти устройства. Ученые издавали статьи, книги, посвященные этой проблеме. Однако существовало непреодолимое препятствие – нестабильность электронных элементов, на которых осуществлялась реализация датчиков случайных чисел. Из-за нестабильности эти датчики необходимо было постоянно настраивать.

Появление высокопроизводительной вычислительной техники и развитых языков программирования упростило задачу создания эталонов. Датчики или генераторы случайных чисел стали создавать программным путем. В библиотеках языков программирования присутствуют подпрограммы (процедуры) с названиями random, randomizir, randu и прочее, которые позволяют получать на ЭВМ псевдослучайные числа, квазиравновероятно распределенные в интервале [0, 1].

Строго говоря, для формирования случайных событий с различными функциями распределения в ЭВМ необходимо иметь, как эталон, совокупность случайных чисел с равномерным законом распределения. Эта совокупность случайных чисел будет эталоном, т.е. базовой последовательностью случайных чисел. На рис. 4.1 показан вид функции распределения случайных чисел, равномерно (равновероятно) распределенных в интервале [0, 1].

Рис. 4.1

Функция распределения вероятностей случайных чисел, равномерно распределенных в интервале [0, 1], имеет вид

      (4.1)

Математическое ожидание (генеральное среднее) случайных чисел, равномерно распределенных в интервале [0, 1], будет равно M[х]=0,5, а дисперсия D[х]=1/12.

Однако применение возможностей ЭВМ не позволяет теоретически получить последовательность случайных чисел с равномерным распределением. Действительно, любое число в ЭВМ формируется в регистре. Если число разрядов регистра ЭВМ равно k, то число Х в десятеричной системе счисления будет сформировано согласно формуле

      (4.2)

в которой коэффициенты i{0, 1} составляют двоичный код k-1,k,…,2,1,0 случайного десятеричного числа Х.

Десятеричное число Х будет случайным и будет иметь равномерное распределение, если коэффициенты i=0 с вероятностью рi=0,5 или i=1 с вероятностью рi=1-0,5. Тогда числа Х={i/(2k-1)} принимают значения i/(2k-1) (i=0,1,2,...,2k-1) с постоянной вероятностью Р=1/2k.

Но полученные при этих условиях случайные числа Х будут иметь не равномерное, а квазиравномерное равновероятностное распределение в интервале [0,1], вид которого показан на рис. 4.2, что связано с дискретностью чисел Х. Математическое ожидание и дисперсия для данного распределения определяются соотношениями:

    (4.3)

.  (4.4)

Рис. 4.2

Таким образом, математическое ожидание квазиравномерного равновероятностного распределения в интервале [0,1] M[х] точно совпадает с генеральным средним для равномерного распределения в интервале [0,1], а дисперсия при k асимптотически стремится к дисперсии для равномерного распределения, равной 1/12. Практически при k>15 обеспечивается требуемая точность в имитационных исследованиях.

В ЭВМ нет генератора случайных чисел i, принимающих значения либо 0, либо 1 с вероятностью Р=1/2, поэтому аппаратная реализация датчика не имеет смысла. Датчики случайных чисел разрабатывают в виде программ, поэтому случайные числа, получаемые на ЭВМ, не являются случайными. Эти числа формируются на основе детерминированных преобразований в виде алгоритмов, поэтому их называют псевдослучайными. Так как алгоритм представляет собой строгую последовательность действий, то последовательность, получаемых псевдослучайных, квазиравномерно распределенных в интервале [0,1] чисел, имеет период Р, например: Х0, Х1, Х2, Х3,...,ХР-2, ХР-1, ХР, Х0, Х1, Х2, Х3,...

Датчики псевдослучайных, квазиравномерно распределенных в интервале [0,1] чисел создают таким образом, чтобы период Р был как можно более большим и превосходил в несколько раз число испытаний (опытов), производимых на ЭВМ в имитационной программе. Если при моделировании число обращений к программному датчику случайных чисел оказывается меньше периода, измеряемого числом различных случайных чисел, то такая периодичность программного датчика не оказывает существенного влияния на результаты моделирования.

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

а) аналитические;

б) методы перемешивания.

При использовании аналитических методов в псевдослучайной последовательности Х0, Х1, Х2, …, Хr-1 очередное число Хr получают с помощью некоторого выбранного рекуррентного соотношения , аргументами которого являются одно или несколько предыдущих чисел последовательности, т.е. Хr=r-1, Хr-2,..., Х0).

Самым простым примером может служить метод вычетов, в котором для получения числа Хi+1 используется следующее рекуррентное соотношение:

Хi+1=bХi(mod M),       (4.5)

где выражение i(mod M) означает остаток от деления произведения i на число M; Хi+1 - очередное случайное число; Хi - предыдущее случайное число; b - некоторая константа; M - число, определяющее значение получаемых случайных чисел.

Рассмотрим пример. Пусть очередное число Хi+1 определяется по формуле

,      (4.6)

где А и В заданные константы, i+1[ - операция взятия мантиссы числа Хi+1. На рис. 4.3 приведен алгоритм генератора псевдослучайных, квазиравномерно распределенных чисел в интервале [0,1].

Рис. 4.3

В блоке 1 алгоритма осуществляется задание начального такта моделирования Т=0, заданное число тактов моделирования (генерации) TZ, равное количеству чисел Х, которое должно быть получено от датчика. В блоках 2-6 вычисляется очередное число Х[T+1]. В блоке 8 проверяется условие генерации датчиком заданного числа чисел Х. Если условие ТTZ выполняется, то наращивается в блоке 9 такт моделирования. Таким образом, генератор псевдослучайных, квазиравномерно распределенных чисел может быть реализован согласно заданной формуле (4.6).

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

На рис. 4.4 показан пример метода перемешивания.

Число Хi вначале циклически сдвигается на три разряда влево (K=3), а затем полученное после этого сдвига число поразрядно суммируется по модулю два с начальным числом Хi. Получаем двоичное число Хi*.

Затем это число циклически сдвигается на два разряда вправо (L=2) и полученный результат суммируется поразрядно по модулю два с числом Хi*. В результате получим последующее число псевдослучайной последовательности Хi+1.

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

Рис. 4.4

Для того чтобы числа были распределены в интервале [0,1], необходимо провести соответствующее нормирование, например, считать, что получаемые числа это дробная часть.

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

В блоке 1 алгоритма осуществляется ввод начального такта моделирования Т=0, заданное число тактов моделирования (генерации) TZ, начального массива Х0[I] размерностью IM, значение IM, а также K - число сдвигов влево и L - число сдвигов вправо. В массиве Х0[I] записана начальная константа Х0.  В блоках 3-6 определяется массив ХK[I], полученный из массива Х0[I] в результате циклического сдвига его элементов на K разрядов влево.

Например, если IM=7, а K=3, то при I=0 (блок 4) ХK[3]=Х0[0] (блок 5), при I=1 ХK[4]=Х0[1], при I=2 ХK[5]=Х0[2], при I=3 ХK[6]=Х0[3], при I=4 ХK[0]=Х0[4], при I=5 ХK[1]=Х0[5], при I=6 ХK[2]=Х0[6]. Результатом операции (I+K)modIM будет остаток от деления числа I+K на число IM.

Рис. 4.5

Затем в блоках 7-10 алгоритма происходит поразрядное суммирование по модулю два элементов массивов ХK[I] и Х0[I]. В результате будет получен массив ХМ[I].

В блоках 11-14 определяется массив ХL[I], полученный из массива ХM[I] в результате циклического сдвига его элементов на L разрядов вправо. Например, при L=2, при I=0 (блок 12) ХL[0]=ХM[2] (блок 13), при I=1 ХL[1]=ХM[3], при I=2 ХL[2]=ХM[4], при I=3 ХL[3]=ХM[5], при I=4 ХL[4]=ХM[6], при I=5 ХL[5]=ХM[0], при I=6 ХL[6]=ХM[1].

Затем в блоках 15-18 алгоритма происходит поразрядное суммирование по модулю два элементов массивов ХМ[I] и ХL[I]. В результате будет получен массив ХT[I], элементы которого представляют собой двоичный код мантиссы искомого псевдослучайного, квазиравномерно распреде-ленного в интервале [0,1] числа. В блоке 19 элементы массива ХT[I] присваиваются элементам массива Х0[I]. В блоке 21 выполняется условие генерации датчиком заданного числа чисел Х.

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

4.3. Проверочные тесты

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

Тест частот. Отрезок [0,1] разбивается на m (обычно 10-20) равных интервалов, как это показано на рис. 4.6.

Рис. 4.6

Датчик псевдослучайных, квазиравномерно распределенных чисел генерирует N величин, каждая из которых принадлежит одному из m отрезков.

В результате проведения N испытаний будут получены эмпирические частоты ni, представляющие собой число попаданий чисел датчика в интервал i, i=. Деление частот ni на число опытов N даст эмпирические частоты.

Если рассматривать теоретическое равномерное распределение случайной величины на отрезке [0,1], то теоретические вероятности попадания в каждый из m отрезков одинаковы и равны значению 1/m.

Для проверки согласия датчика псевдослучайных, квазиравномерно распределенных чисел и теоретического равномерного распределения полученные эмпирические частости ni/N, (i=),  сравнивают с теоретическими вероятностями 1/m. Согласие проверяется по критерию 2, т.к. случайная величина

     (4.7)

подчиняется распределению 2 с (m-1) степенями свободы, где N - объем выборки (число опытов).

На рис. 4.7 приведен алгоритм теста частот. Рассмотрим его работу.

Рис. 4.7

В подпрограмме WWOD осуществляется ввод исходных данных для моделирования: N=0 - начальный такт моделирования; NZ - заданное число тактов моделирования (число опытов); m – число интервалов разбиения отрезка [0,1]. В подпрограмме WWOD также осуществляется обнуление необходимых для работы идентификаторов и счетчиков.

В блоке 2 происходит наращивание тактов моделирования. В блоке 3 датчиком случайных чисел генерируется число Х. Затем это число Х сравнивается с правыми границами интервалов разбиения отрезка [0,1]. Для этого в блоках 4 - 6 организован цикл по переменной I и сравнение числа Х с числом А, которое последовательно принимает значения: 1/m, 2/m, 3/m, …, m/m.

При выполнении условия ХA содержимое соответствующего счетчика увеличивается на единицу (см. блок 7).

В блоке 8 осуществляются проверки выполнения числа опытов. Если датчик псевдослучайных, квазиравномерно распределенных чисел выполнил генерацию заданного числа опытов NZ, то в блоке 9 происходит вычисление случайной величины 2.

В подпрограмме WIWOD (см. блок 10) на экран дисплея выводятся значения случайной величины 2, счетчиков K[I], в которых подсчитаны частоты ni, а также могут быть выведены гистограммы эмпирических частостей =ni/N=K[I]/NZ, как это показано на рис. 4.8.

На рис. 4.8 приведены гистограммы для заданного числа опытов NZ =10000. Очевидно, что форма гистограммы частот должна совпадать с формой гистограммы частостей, т.к. эмпирические частости определяются по формуле =K[I]/NZ.

Рис. 4.8

Тесты пар частот. Пусть датчик псевдослучайных, квазиравномерно распределенных чисел генерирует последовательность чисел Х12,...,ХN. Рассматриваются последовательные пары случайных чисел. Квадрат со сторонами [0,1] на [0,1] делится на m2 частей, как это показано на рис. 4.9.

Если пары образовать в виде 12), (Х34)..., то каждая пара случайно попадает в одно из m2 делений квадратной таблицы. Пары 12), (Х34)... взаимно независимы и результат их попадания в одно из m2 делений квадратной таблицы оценивается эмпирической частотой nij.

Рис. 4.9

Если рассматривать теоретическое равномерное распределение случайной величины на отрезке [0,1] и образование из чисел таких же пар, то теоретические вероятности попадания в каждый из m2 делений квадратной таблицы одинаковы и равны значению 1/m2.

Для проверки согласия по данному тесту пар частот датчика псевдослучайных, квазиравномерно распределенных чисел и теоретического равномерного распределения полученные эмпирические частости 2nij/N, (i,j=),  сравнивают с теоретическими вероятностями 1/m2. Согласие проверяется по критерию 2, т.к. случайная величина

  (4.8)

распределена по закону 2 с (m2-m) степенями свободы, где N/2 - объем выборки пар случайных величин, генерированных датчиком псевдослучайных, квазиравномерно распределенных чисел. На рис. 4.10 приведен алгоритм для рассмотренного теста пар частот.

Рис. 4.10

В подпрограмме WWOD осуществляется ввод исходных данных для моделирования:

N=0 - начальный такт моделирования;

NZ - заданное число тактов моделирования;

m - число интервалов разбиения отрезка [0,1], происходит обнуление необходимых для работы идентификаторов и счетчиков.

В блоке 2 происходит наращивание тактов моделирования. В блоке 3 датчиком случайных чисел генерируется число Х. Затем в блоках 4 – 6 определяется индекс I. В блоке 7 датчиком случайных чисел генерируется второе число Х (тем самым образована пара 12)). В блоках 8 – 10 определяется индекс J.

В блоке 11 в счетчиках K[I,J] осуществляется подсчет частот попадания случайных пар kр) в соответствующие деления квадратной таблицы. В блоке 12 осуществляются проверки выполнения числа опытов. Если датчик псевдослучайных, квазиравномерно распределенных чисел выполнил генерацию заданного числа опытов NZ/2, то в блоке 13 происходит вычисление случайной величины 2. В подпрограмме WIWOD (см. блок 14) на экран дисплея выводятся значения случайной величины 2, счетчиков K[I,J], значения эмпирических частостей =2nij/N=2K[I,J]/NZ. Генерируемая выборка в данном методе используется неэффективно. Можно пары образовать в следующем виде 12),(Х23),(Х34),... .

Этот метод образования пар более эффективен, т.к. полнее использует выборку чисел, но из-за зависимостей пар случайная величина 2 определится по формуле

],   (4.9)

где .

Случайная величина 2 будет определена распределением 2 с (m2-m) степенями свободы.

На рис. 4.11 приведен алгоритм для рассмотренного теста пар частот с эффективным использованием выборки.

Рис. 4.11

В данном алгоритме в блоке 2 генрируется число Х1, а затем в блоках 3 - 5 определяется индекс I. В блоке 6 происходит наращивание тактов моделирования. В блоке 7 датчиком случайных чисел генерируется число Х. В блоках 8 – 10 определяется индекс J (тем самым образована пара 12)). В блоке 11 в счетчиках K[I,J] осуществляется подсчет частот попадания случайных пар в деления квадратной таблицы. В блоке 12 индексу I присваивается значение индекса J. При N=2 будет образована пара 23) и т.д.

4.4. Имитация случайных событий

Пусть события S1, S2,..., Sm образуют полную группу несовместимых событий, каждое из которых может произойти с вероятностью Рi, причем . Это означает, что в результате одного опыта может произойти лишь одно из событий S1, S2,..., Sm. Рассмотрим метод и алгоритм имитации данных случайных событий.

Разбиваем отрезок [0,1] на m частей длиной Р1,Р2,...,Рm, при этом точки деления отрезка имеют следующие координаты:

,

как это показано на рис. 4.12.

Рис. 4.12

Имитация осуществляется следующим образом. Пусть теперь Х - очередное число от датчика псевдослучайных, квазиравномерно распределенных чисел, который далее для упрощения будем называть генератором случайных чисел. Если число генератора Аk-1Х<Аk, то считаем, что произошло событие Sk. Действительно

Р(Sk)=Р(Ak-1<Х<Ak)=Ak-Ak-1=Рk.

Алгоритм имитации несовместимых событий S1, S2,..., Sm, образующих полную группу, рассмотрим для случая m=5. Схема алгоритма приведена на рис. 4.13.

Рис. 4.13

На каждом такте моделирования случайное число Х сравнивается со значениями А1, А2, А3, А4, А5 в блоках 4, 6, 8, 10 соответственно.

Если случайное число ХА1, то произошло событие S1 и содержимое счетчика K[1] увеличивается на единицу. Если случайное число ХА2, то произошло событие S2 и содержимое счетчика K[2] увеличивается на единицу и т.д. События S1, S2,..., Sm могут образовывать неполную группу несовместимых событий. Если событие Si происходит с вероятностью Рi, то для этих событий . 

Для имитации событий S1, S2,..., Sm, образовывающих неполную группу, необходимо ввести дополнительное событие S0, с которым связано непоявление ни одного из вышеперечисленных событий. Вероятность события S0 определится по формуле . После введения события S0 имитация неполной группы событий осуществляется по аналогии с рассмотренной выше имитацией событий, образующих полную группу несовместимых событий.

Рассмотрим ситуацию, когда осуществляется имитация зависимых событий. Событие В будет зависеть от события А, если после наступления события А событие В может произойти с условной вероятность Р(В/А).

Метод и алгоритм имитации рассмотрим на примере.

Пусть события S1, S2, S3, S4 образуют полную группу независимых событий. Если происходит событие S1, то с вероятностью Р(U1/S1)=Р(11) может произойти событие U1. Если происходит событие S2, то с вероятностью Р(H1/S2)=Р(12) может произойти событие H1 и с вероятностью Р(H2/S2)=Р(22) может произойти событие H2, причем события H1 и H2 образуют полную группу событий. Если происходит событие S4, то с вероятностью Р(D1/S4)=Р(14) может произойти событие D1.

На рис. 4.14 приведен граф событий.

Рис. 4.14

На рис. 4.15 приведена схема алгоритма имитации рассмотренной системы случайных событий.

После генерации случайного числа Х (см. блок 3) в блоке 4 проверяется условие появления события А1.

Если событие А1 произошло, то содержимое счетчика KS[1] увеличивается на единицу (см. блок 5) и генерируется случайное число Х (см. блок 6). В блоке 7 проверяется условие появления события U1. Если событие U1 произошло, то содержимое счетчика KU[1] увеличивается на единицу (см. блок 8). Если событие А1 не произошло, то в алгоритме выполняется переход от блока 4 к блоку 9 и проверяется условие появления события А2. Если событие А2 произошло, то содержимое счетчика KS[2] увеличивается на единицу (см. блок 10) и генерируется случайное число Х (см. блок 11).

В блоке 12 проверяется условие появления события H1. Если событие H1 произошло, то содержимое счетчика KH[1] увеличивается на единицу (см. блок 13). Если событие H1 не произошло, то считается, что произошло событие H2, т.к. события H1 и H2 составляют полную группу. Содержимое счетчика KH[2] увеличивается на единицу (см. блок 14).

Рис. 4.15

Если событие А2 не произошло, то в алгоритме выполняется переход от блока 9 к блоку 15 и проверяется условие появления события А3. Если событие А3 произошло, то содержимое счетчика KS[3] увеличивается на единицу (см. блок 16).

Если событие А3 не произошло, то считается, что произошло событие А4. Содержимое счетчика KS[4] увеличивается на единицу (см. блок 17) и генерируется случайное число Х (см. блок 18). В блоке 19 проверяется условие появления события D1. Если событие D1 произошло, то содержимое счетчика KD[1] увеличивается на единицу (см. блок 20). В блоке 21 осуществляется контроль выполнения цикла по переменной N.

4.5. Имитация непрерывных случайных величин

Если событие Х принимает значения в некоторой области непрерывных величин, то для аналитического моделирования непрерывных событий применяют функцию распределения вероятностей F(Х<х) или плотность распределения вероятностей f(х). Функция распределения вероятностей F(Х<х) определяет вероятность того, что событие (случайная величина) Х меньше либо равна некоторому значению х, т.е.

F(Х<х)=Р{Х<х}.

Вид функции распределения вероятностей приведен на рис. 4.16. Это монотонно возрастающая (неубывающая) функция в диапазоне от нуля до единицы, не имеющая разрывов на интервале значений от + до +.

Плотность распределения вероятностей f(х), называемая еще дифференцированным распределением вероятностей, определяется по формуле

, причем .

Рис. 4.16

Гипотетический вид плотности распределения вероятностей f(х) приведен на рис. 4.17.

Рис. 4.17

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

- экспоненциальное распределение

F(х)=1-e-х;

- распределение Пуассона, определяющее вероятность появления k событий за время t по формуле

,

где - математическое ожидание;

- распределение Эрланга r-го порядка, плотность распределения вероятностей для которого определится по формуле

;

- нормальное (гауссово) распределение, плотность распределения вероятностей для которого определится по формуле

,

где m – математическое ожидание, - среднеквадратичное отклонение, а также другие распределения.

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

4.5.1. Метод обратной функции. Пусть случайная величина Х определена функцией распределения вероятностей F(х) и плотностью распределения вероятностей f(х).

Если Р - случайная величина, равномерно распределенная на отрезке [0,1], то случайная величина Х может быть получена из решения следующего уравнения:

.      (4.10)

На рис. 4.18 приведена иллюстрация метода обратных функций. Как следует из рис. 4.18, выполняется как бы обратное преобразование вероятности Р в случайную величину Х, поэтому данный метод и получил наименование метода обратных функций.

Рис. 4.18

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

- вырабатывается датчиком случайной равномерной последовательности случайное число Р;

- случайная величина (случайное число) Х, имеющее распределение f(х), находится из решения уравнения (4.10).

Рассмотрим пример.

Пусть случайная величина Х определена функцией распределения вероятностей F(х)=1-e-х. Плотность распределения вероятностей f(х)= e-х. При известной вероятности Р значение случайной величины Х определится по формуле

.

Таким образом, получили трансцендентное уравнение Р=1-e-Х с одним неизвестным Х. Решение этого уравнения:

.

В силу того, что Р – число, равномерно распределенное на отрезке [0,1], то и 1-Р – число, равномерно распределенное на отрезке [0,1], поэтому случайную величину Х находят по формуле

.

На рис. 4.19 приведена схема алгоритма метода обратных функций для рассмотренного случая.

Рис. 4.19

Подпрограммы WWOD и WIWOD предназначены для реализации интерфейса пользователя и инициализации программного модуля. Блоками 2 и 6 организован цикл по переменной N. В подпрограмме GEN генерируется число Р, равномерно распределенное на отрезке [0,1]. В блоке 4 определяется случайная величина Х, определенная функцией распределения вероятностей F(х)=1-e-х. Подпрограмма STAT предназначена для набора и обработки статистических данных о значениях случайной величины Х.

4.5.2. Метод ступенчатой аппроксимации. Метод обратных функций применим для моделирования непрерывных случайных величин Х в том случае, если существует функция распределения вероятностей F(х), т.е. уравнение 10 имеет аналитическое решение. Однако известны распределения случайных величин, например, нормальное распределение, для которых функция распределения вероятностей F(х) аналитически не определяется.

Для имитации случайных величин Х, определяемых только плотностью распределения вероятностей f(х) применяется метод ступенчатой аппроксимации. Рассмотрим суть этого метода. Зависимость плотности распределения f(х) представляется графически в интервале изменения случайной величины Х от a до b.

Если случайная величина задана на [-,+], то достаточно ограничиться минимально и максимально возможными числами, воспринимаемыми ЭВМ.

Суть метода ступенчатой аппроксимации показана на рис. 4.20.

Разобьем [a,b] на n интервалов таким образом, чтобы площади под кривой f(х) внутри каждого интервала были равны, как это показано на рис. 4.20, т.е.

,

где Сi () - координаты точек разбиения.

Рис. 4.20

Вероятность того, что случайная величина Х попадет в любой из интервалов i-1, Сi], , определится по формуле

,

т.е. попадание на любой отрезок равновероятно.

Если внутри интервала распределение случайной величины Х также равномерное, то значение случайной величины Х на оси х может быть определено, как Х=Сi-1+, где  - равномерно распределенная случайная величина на интервале i-1, Сi], представляющая собой расстояние от левого конца (Сi-1) i–го интервала. На рис. 4.21 приведена схема алгоритма генерации случайных величин Х с применением метода ступенчатой аппроксимации.

Рис. 4.21

Правило имитации случайных величин Y сводится к следующему:

- получаем от генератора равномерно распределенных чисел случайное число Р1 (см. блок 3);

- из значения числа Р1 находим индекс i=]nР1[ интервала i-1, Сi], где ]nР1[ - целая часть числа 1, причем ]nР1[<nР1 (см. блок 4);

- получаем от генератора равномерно распределенных чисел случайное число Р2 (см. блок 5);

- из значения числа Р2 находим случайную величину 2ii-1), т.е. значение числа Р2 «приводится» к величине интервала i-1, Сi] (см. блок 6;

- находим случайную величину Х, имеющую заданную плотность распределения вероятностей f(х), по формуле (см. блок 7)

.

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

4.5.3. Использование предельных теорем. Для имитации случайной величины Х, имеющей нормальный закон распределения вероятностей, используют свойство сходимости независимых величин к нормальному распределению. Для получения нормального распределения чисел с параметрами: математическое ожидание mХ=0, среднеквадратичное отклонение Х=1 удобен искусственный прием, основанный на центральной предельной теореме теории вероятностей. На рис. 4.22 приведен алгоритм получения случайный величин Х с применением свойств центральных предельных теорем.

Согласно центральной предельной теореме при достаточно большом значении n величина Z может считаться нормально распределенной с параметрами

, .

Для имитации случайной величины Х в качестве исходных чисел возьмем k равномерно распределенных на отрезке [1,-1] случайных чисел, получаемых из интервала [0,1] по правилу zi=2Рi-1 (см. блоки 5,6).

Рис. 4.22

Сформируем величину Z согласно следующей формуле (см. блоки 3 – 7):

.

Выполним нормирование величины Z и получим (см. блок 9)

.    (4.11)

Случайная величина Х будет иметь нормальное распределение с mХ=0, Х=1.

Установлено, что при k>8 формула (4.11) дает хорошие результаты.

4.6. Имитация марковского процесса

4.6.1. Моделирование дискретной цепи Маркова. Рассмотрим дискретную цепь Маркова или марковский процесс с дискретным временем перехода из одного состояния в другое. Математическая модель дискретной цепи Маркова задается матрицей вероятностей переходов ||Рij||, имеющей вид

,     (4.12)

где Рij - вероятность перехода из состояния zi в состояние zj в некоторый дискретный момент времени kt (k=1,2,3,...).

Начальное состояние марковского процесса определяется матрицей-строкой начальных вероятностей ||Р0||=|Р1(0), Р2(0), ..., Рn(0)|, где Рi(0) -  вероятность нахождения процесса в zi-м состоянии при t=0. Вероятности перехода Рij не зависят от времени, т.е. процесс является однородным.

На рис. 4.23 приведен обобщенный алгоритм имитации дискретной цепи Маркова.

Рис. 4.23

Подпрограммы WWOD и WIWOD реализуют интерфейсную часть имитационной программы. В подпрограмме WWOD осуществляется ввод элементов массива Р0[I], в который заносятся вероятности Рi(0), и массива Р[I,J], в который заносятся вероятности Рij. Определяется начальный такт моделирования T=0 и заданное число тактов моделирования TZ.

Подпрограмма OРRZ0 предназначена для определения начального состояния, а подпрограмма OРRZ – для определения состояний в процессе моделирования смены состояний. Подпрограмма STAT предназначена для набора статистических данных.

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

Определяем числовые границы

.

Исходя из значения случайного числа Р0, генерированного датчиком случайных чисел, определяется номер r начального состояния zr(0), для которого будет справедливо условие

.

Затем датчик случайных чисел вырабатывает случайное число Р1, которое также сравнивается с границами

.

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

Рассмотрим реализации подпрограмм.

На рис. 4.24 приведен алгоритм подпрограммы OРRZ0.

На рис. 4.25 приведен алгоритм подпрограммы OРRZ.

На рис. 4.26 приведен алгоритм подпрограммы STAT.

В алгоритме подпрограммы OРRZ0 (см. рис. 4.24) в блоке 1 вырабатывается число Р датчиком случайных чисел.

Затем реализуется цикл по переменной I для сравнения числа Р с элементами массива Р0[I]. Для этого введен идентификатор А, который в блоке 2 определен А=0.

Рис. 4.24

Рис. 4.25

Рис. 4.26

При наращивании переменной I идентификатор А будет принимать последовательно значения: I=1, А=Р1(0); I=2, А=Р1(0)+Р2(0); I=3, А=Р1(0)+Р2(0)+Р3(0); …; I=n, А=Р1(0)+Р2(0)+Р3(0)+…+Рn(0).

При первом же выполнении условия РA считается, что найден индекс начального состояния zi(0). Таким образом, выходным параметром подпрограммы OРRZ0 является значение индекса I, при котором выполнено условие РА=Р1(0)+Р2(0)+Р3(0)+…+Рi(0).

В алгоритме подпрограммы OРRZ (см. рис. 4.25) в блоке 1 датчиком случайных чисел формируется случайное число Р[0,1]. Затем реализуется цикл по переменной J для сравнения числа Р с элементами массива Р[I,J]. Введен идентификатор В. При наращивании переменной J идентификатор В будет принимать последовательно значения: J=1, А=Р(I,1); J=2, А=Р(I,1)+Р(I,2); …; J=n, А=Р(I,1)+Р(I,2)+Р(I,3)+…+Р(I,n).

При выполнении условия РB считается, что найдет индекс состояния zj(T) в текущем такте моделирования T.

Выходным параметром подпрограммы OРRZ является значение индекса J, при котором выполнено условие РВ=Р(I,1)+Р(I,2)+Р(I,3)+…+Р(I,J). 

В алгоритме подпрограммы STAT (см. рис. 4.26) в блоке 1 в счетчиках K[J] осуществляется подсчет частот появления событий zj. Затем в блоке 2 осуществляется присвоение значения индекса J индексу I, т.е. на следующем такте моделирования в подпрограмме OРRZ будет выполнен анализ j-й сроки матрицы вероятностей переходов ||Рij||. Алгоритм имитации дискретной цепи Маркова может быть представлен в более сокращенном виде, как показано на рис. 4.27. Объединим матрицы ||Р0|| и ||Рij|| в обобщенную матрицу вида

.     (4.13)

В подпрограмме WWOD осуществляется ввод элементов массива Р[I,J], в который заносятся вероятности Рij, причем , . Определяется начальное значение индекса I=0, начальный такт моделирования T=0 и заданное число тактов моделирования TZ.

Так как в первом такте моделирования (Т=1) I=0, то в первом такте будет рассматриваться верхняя строка обобщенной матрицы ||Рij||, т.е. фактически матрицы ||Р0||.

В блоке 3 формируется случайное число Р[0,1]. Затем реализуется цикл по переменной J для сравнения числа Р с элементами массива Р[I,J] при I=0 (см. блоки 4 – 6). При выполнении условия РВ в соответствующий счетчик K[J] добавляется единица (см. блок 7). В блоке 8 индексу I присваиваются значения индекса J. Процесс моделирования продолжается до тех пор, пока не будет выполнено условие окончания моделирования T<TZ.

Рис. 4.27

В счетчиках K[J] после окончания моделирования будут подсчитаны частоты пребывания марковского процесса в состояниях zj, . В подпрограмме WIWOD осуществляется вывод значений счетчиков K[J], а также вывод эмпирических оценок финальных вероятностей марковского процесса, определяемых по формуле

.

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

- распределением вероятностей начального состояния процесса Рi(0) в момент t0;

- матрицей вероятностей переходов ||Рij|| (12);

- матрицей-строкой функций распределения времени пребывания в состоянии |A|={A1(t), A2(t), …, An(t)}, где Ai(t) - функция распределения времени пребывания марковского процесса в состоянии zi.

На рис. 4.28 приведен обобщенный алгоритм имитации вложенной цепи Маркова.

Рис. 4.28

Подпрограммы WWOD и WIWOD предназначены для реализации интерфейса имитационной программы. В подпрограмме WWOD осуществляется ввод элементов массива Р0[I], массива Р[I,J], а также параметров функций распределения Ai(t). Определяется T=0 и TZ. В подпрограмме OРRZ0 (см. блок 3) определяется начальное состояние, в подпрограмме OРRZ определяются состояния в процессе моделирования смены состояний.

Подпрограмма STAT1 предназначена для набора статистических данных о частотах пребывания марковского процесса в состояниях.

В подпрограмме OРRTAU определяется время пребывания марковского процесса в состоянии, определенном в подпрограмме OРRZ. Подпрограмма STAT2 предназначена для набора статистических данных о времени пребывания марковского процесса в состояниях.

Рассмотрим реализацию подпрограмм OРRZ0, OРRZ, OРRTAU.

Подпрограммы OРRZ0, OРRZ имеют такой же вид, как и аналогичные подпрограммы для дискретной цепи Маркова (см. рис. 4.24 и рис. 4.25). Выходным параметром подпрограммы OРRZ является индекс J состояния zj(T) в текущем такте моделирования T.

Схема алгоритма подпрограммы OРRTAU будет иметь вид алгоритмов генерации случайной величины, рассмотренных в разд. 4.5, метод обратных функций, метод ступенчатой аппроксимации, использование предельных теорем. О схемной реализации подпрограмм STAT1 и STAT2 будет сказано ниже.

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

- распределением вероятностей состояния Рi(0) в момент t0;

- матрице , имеющей вид

,    (4.14)

где ij - интенсивностей переходов марковского процесса из состояния zi в состояние zj.

Способ моделирования сводится к следующему.

Пусть в момент t состояние марковского процесса z(t)=zi. Известен момент времени ti, в который марковский процесс должен покинуть состояние zi. Для определения следующего состояния zi+1 и момента ti+1 перехода в последующее состояние генерируется ряд независимых случайных величин ij по показательному закону распределения с параметрами i-й строки матрицы ||ij ||. Затем определяется минимальное значение из полученной совокупности независимых случайных величин ij, т.е.

.     (4.15)

Следующий момент изменения состояния будет ti+1=ti+ij, а состояние, в которое переходит система в момент ti+1, будет zj, где zj - состояние, при котором zj было минимальным из всех значений совокупности независимых случайных величин {ii}. Физическая интерпретация такого подхода к моделированию марковского процесса с непрерывным временем позволяет разрабатывать имитационные модели реальных систем.

Схема алгоритма моделирования вложенной цепи Маркова по данному способу приведена на рис. 4.29.

Рис. 4.29

Отличие алгоритма по схеме развернутой рекуррентной имитации (см. рис. 4.29) от алгоритма, приведенного на рис. 4.28, состоит в блоках 4 и 5, в которых реализована генерация случайных величин ki и формулы (4.15).

На рис. 4.30 приведена схема алгоритма подпрограммы OРRTAU для данного способа. На рис. 4.31 приведена схема алгоритма подпрограммы MINTAU.

Рис. 4.30

В подпрограмме OРRTAU (см. рис. 4.30) входным параметром является индекс I предыдущего состояния марковского процесса.

Для каждого значения J в блоке 4 по методу обратных функций определяется значение случайной величины TAU[J], отождествляющей собой описанную выше величину ii. Таким образом, в данной подпрограмме будет сформирован массив TAU[J].

В подпрограмме MINTAU (см. рис. 31) входным параметром являются элементы массива TAU[J], из которых необходимо выбрать минимальный элемент.

В блоках 1 – 5 осуществляется выбор минимальной величины (MIN) из массива TAU[J]. Происходит это следующим образом. При переменной массива J=1 первый элемент массива TAU[1] определяется, как минимальный элемент MIN (см. блок 1).

Рис. 4.31

Затем второй элемент массива TAU[2] сравнивается со значением MIN (см. блок 3). Если второй элемент массива TAU[2] меньше значения MIN, то идентификатору MIN присваивается значение TAU[2] (см. блок 4). Если условие TAU[2]<MIN не выполняется, то рассматривается условие TAU[3]<MIN для элемента массива TAU[3] и т.д.

В блоках 6 – 8 определяется индекс J элемента массива TAU[J], равного числу MIN. Это индекс J определяет индекс состояния zJ, в который переходит марковский процесс.


5. ОБРАБОТКА РЕЗУЛЬТАТОВ МОДЕЛИРОВАНИЯ НА ЭВМ

5.1. Выбор числа опытов

При разработке имитационных моделей для исследования случайных объектов существует задача выбора числа опытов (объема выборки). Это непростая задача, т.к. во-первых, необходимо обосновать достоверность результатов моделирования и связать достоверность с точностью, а во-вторых, существуют события, вероятность появления которых очень мала (Р0) или, наоборот очень велика (Р1). Для обоснования объема выборки при имитационном моделировании применяют аналитические подходы [8,12,13].

В инженерной практике известны критерии для оценки погрешности. Непрерывную (аналоговую) величину х(t) в ряде инженерных задач рассматривают как дискретную, как показано на рис. 5.1.

Рис. 5.1

Если х*(t) - результат измерения непрерывной величины х(t), то для любого момента t текущая погрешность дискретизации этой непрерывной величины определится: (t)=х(t)-х*(t). Выбор критерия оценки (t) зависит от назначения величины х(t). Известны следующие критерии.

Критерий наибольшего отклонения имеет вид

. Критерий применим, если известны априорные сведения о сигнале в форме условия Липшица , где l -некоторая константа.

Среднеквадратичный критерий приближения определяется по формуле

.

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

Интегральный критерий как мера отклонения х(t) от х*(t) имеет вид

.

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

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

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

|A-х*|<,        (5.1)

где - точность оценки.

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

Р(|A-х*|<e)=.       (5.2)

Воспользуемся сформулированным критерием (5.2) для определения точности результатов методом статистического моделирования. Пусть цель моделирования - вычисление вероятности р появления события A. Количество наступления события A в реализации процесса является случайной величиной, принимающей значение х1=1 с вероятностью р и значение х2=0 с вероятностью 1-р. Пример случайного появления события A показан на рис. 5.2.

Рис. 5.2

Математическое ожидание случайной величины определится по формуле

M[]=х1р+х2(1-р)=р.      (5.3)

Если х1=1, как появление события A, а х2=0, как не появление события, то значение M[] совпадает с вероятностью р наступления события A.

Дисперсия определится по формуле

D[]=[х1-M[]]2р+[х2-M[]]2(1-р)=р(1-р).  (5.4)

При выполнении имитационного моделирования оценкой вероятности р является частость m/N наступления события A при N реализациях, m - число испытаний, в которых событие A наступило. Частость m/N определяется формулой

,    (5.5)

где i - количество наступлений событий A в реализации с номером i.

Из формул (5.3), (5.4) и (5.5) можно определить математическое ожидание и дисперсию частости m/N 

.      (5.6)

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

       (5.7)

Подставив в формулу (5.7) значение D из формулы (5.6), получим

       (5.8)

Из формулы (5.8) можно определить количество реализаций N, необходимых для получения оценки m/N с точностью и достоверностью :

     (5.9)

В формуле (5.9) неизвестны величины N и р, т.к. вероятность р определяется, исходя из ее оценки m/N, а N - число необходимых для этого опытов. Поэтому в практике моделирования для определения N поступают следующим образом. Выбирают N0=50-100. По результатам N0 реализаций определяют m/N0, т.е. осуществляют примерную оценку вероятности р=m/N0. Затем по формуле (5.9) окончательно выбирают N, принимая р=m/N.

Другим случаем является оценка по результатам моделирования среднего значения некоторой случайной величины. Пусть непрерывная случайная величина A имеет среднее значение  и дисперсию 2. В реализации с номером i случайная величина A принимает значение хi. В качестве оценки для среднего значения (математического ожидания) A используется среднее арифметическое

.

В силу центральной предельной теоремы при N   будет иметь приблизительно нормальное распределение с математическим ожиданием  и дисперсией 2/N, поэтому точность определится по формуле

.

Число реализаций определится по формуле

.      (5.10)

Так как в формуле (5.10) неизвестными являются число реализаций N и среднеквадратичное отклонение 2, то также выбирают N0=50-100. По результатам N0 реализаций определяют оценку дисперсии, а затем по формуле (5.10) окончательно выбирают N. Количество реализаций N в формуле (5.9) зависит от р, а в формуле (5.10) от 2.

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

5.2. Значимость оценки

5.2.1. Статистическая проверка гипотезы относительно вероятности. При имитационном моделировании получают практические (эмпирические) результаты, которые затем аппроксимируют известными теоретическими распределениями, т.е. выдвигают гипотезу, то данное теоретическое распределение аппроксимирует эмпирическое распределение. Например, необходимо проверить гипотезу относительно того, что при выполнении имитационного моделирования частость р*=m/N является оценкой вероятности р события.

Пусть проверяется настройка станка на среднюю точку поля допуска. Проведено 280 независимых испытаний и интересующее нас событие появилось 151 раз. Модель появления независимых событий – биноминальное распределение. Если гипотетическая вероятность события р=1/2, то математическое ожидание равно =280(1/2)=140, а среднеквадратичное отклонение биноминального распределения определится из формулы =8,37, где q=1-р – вероятность не появления события. Надо получить ответ на вопрос: можно ли считать наблюденную частоту 151 достаточно близкой к теоретической норме 140, отвечающей гипотезе р=1/2.

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

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

На практике, в соответствии с критерием (5.2), задают уровень значимости , т.е. вероятность практически невозможных отклонений. Эта вероятность обычно не превышает значение 0,05. Область больших отклонений, соответствующую уровню значимости , называют критической областью, а само правило проверки – критерием значимости. Критическую границу для отклонений от теоретической нормы можно определить, пользуясь нормальным приближением к биноминальному закону. На рис. 5.3 приведены двухсторонние критические границы для проверки гипотезы р=1/2.

Рис. 5.3

Вероятность 0,95 соответствует при нормированном нормальном распределении интервалу (-1,96, +1,96) около центра распределения, т.к. вероятность того, что абсолютная величина нормированного отклонения превысит значение 1,96, равна 0,05, т.е.

Р(|A|>1,96)=Р(|m-Nр|>1,96)0,05,

где m – практическая частота.

Уровень значимости 0,01 соответствует границе 2,58, т.е.

Р(|A|>2,58)=Р(|m-Nр|>2,58)0,01.

Для рассмотренного выше примера настройки станка на среднюю точку поля допуска =8,37 и 5% критическая граница соответствует 1,968,37=16,41, а 1% критическая граница соответствует 2,588,37=21,59.

Таким образом, область допустимых значений при 5% критической границы определяется пределами

Nр1,96=14016,41,

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

Nр2,58=14021,6.

Если выдвинутая гипотеза, что наблюденная частота 151 достаточно близка к теоретической норме 140, отвечающей вероятности р=1/2, верна, то отклонение частоты от теоретической нормы в пяти случаях из 100 может превышать 16,4, и в одном случае из 100 может превышать 21,6.  Так как в рассмотренном примере отклонение составило 151-140=11, т.е. оно находится в области допустимых значений, то нет оснований считать гипотезу р=1/2 противоречащей наблюдениям.

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

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

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

Для проверки гипотезы согласно критерию выбираются надлежащие уровни значимости (см. разд. 5.2.1) =5%, 2%, 1% и т.д., отвечающие событиям, которые при данном исследовании считаются практически невозможными. Затем определяется критическая область (см. рис. 5.3) данного критерия, вероятность попадания в которую в точности равна уровню значимости , если гипотеза верна. Значения критерия, лежащие вне критической области, образуют дополнительную к ней область допустимых значений (незаштрихованная область на рис. 5.3).

Если /100 – уровень значимости, то вероятность попадания критерия в область допустимых значений при справедливости выдвинутой гипотезы равна 1-/100. Если значение критерия, вычисленное по произведенным наблюдениям (опытам), окажется в критической области, то гипотеза отвергается. Если значение критерия окажется в области допустимых значений, что наблюденное значение критерия не противоречит гипотезе.

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

При заданном уровне значимости можно по разному устанавливать критическую область, гарантирующую этот уровень. Например, в качестве критерия рассматривается некоторый показатель, распределенный при проверяемой гипотезе нормально с плотностью распределения f(x;a;). В качестве критической области соответствующей уровню значимости =5% можно принять:

- область больших положительных отклонений так, что

Р(tq)=Р(x>a+tq)=0,05,

но

,  (5.11)

тогда из таблицы значений нормированной функции Лапласа

определим, что tq=1,65;

- область больших отрицательных отклонений

Р2(tq)=Р(x<a-tq)=Р(x<a-1,65);

 - область больших по абсолютной величине отклонений

Р3(tq)=Р(|x-a|>tq)=0,05,

определив t из соотношения

,  (5.12)

так что tq=1,96;

- область малых по абсолютной величине отклонений

Р4(tq)=Р(|x-a|>tq)=2Ф0(z),    (5.13)

так что tq0,063.

Эти области показаны на рис. 5.4.

1 – область больших положительных отклонений;

2 - область больших отрицательных отклонений;

3 - область больших по абсолютной величине отклонений (состоит из двух половин); 4 - область малых по абсолютной величине отклонений

Рис. 5.4

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

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

Рассмотрим наиболее употребительный критерий 2 (критерий Пирсона). Пусть гипотеза предполагает вид функции распределения Р(х). Вся область изменения случайной величины Х разбита на конечное число k множеств 1, 1, …, k. Если случайная величина Х непрерывна, то множества 1, 1, …, k представляют собой интервалы, а если случайная величина Х дискретна, то множества 1, 1, …, k представляют собой группы отдельных значений случайной величины Х. Пусть pi - вероятность того, что значения случайной величины Х при данном распределении Р(х) принадлежат интервалу i.

Объем выборки N, а mi - число значений случайной величины Х в выборке O(x1, x2, x3, …, xN), попавших в интервал i. Очевидно, что

p1+p2+ …+pk=1,       (5.14)

m1+m2+ …+mk=N.      (5.15)

Если проверяемая гипотеза верна, то mi представляет частоту появления события, имеющего в каждом из N произведенных испытаний вероятность pi. В таком случае mi можно рассматривать как случайную величину, подчиненную биномиальному закону распределения с центром в центре в точке Npi и средним квадратическим . Если N достаточно велико, то можно считать, что частота распределена асимптотически нормально с центром в центре в точке Npi и средним квадратическим .

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

, (i=1,2,…,k),    (5.16)

связанные между собой соотношением

,      (5.17)

вытекающем из условий (5.14) и (5.15).

В качестве меры расхождения данных выборки (эмпирических частот) m1m2, …, mk с теоретическими частотами Np1Np2, …, Npk рассмотрим величину

.     (5.18)

Для практических приложений можно применять подобное равенство:

.       (5.19)

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

Из теории вероятностей известна теорема. Если проверяемая гипотеза верна, то критерий 2, определяемый по формуле (5.18), имеет распределение, стремящееся при N к распределению 2 с k-1 степенями свободы.

При проведении проверки задают уровень значимости % для критерия.

Пусть  обозначает % предел для закона распределения 2 с k-1 степенями свободы. Этот закон имеет табличное задание и его значения приводятся в приложениях книг с изложением теории вероятностей.

Если гипотеза верна, то при достаточно большом числе опытов N справедливо определение вероятности

.      (5.20)

После определения случайной величины 2 по данным выборки O(x1, x2, x3, …, xN), будет выполняться одно из двух условий:

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

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

Во втором случае в % всех случаев, но неизвестно каких, гипотеза неверна. Принято считать достаточным нормальное приближение для практических расчетов, если Npi10 i. Если есть группы со значениями Npi меньшими 10, то рекомендуют соседние группы объединять так, чтобы новые группы удовлетворяли условию Npi10 i.

Если число степеней свободы k>30, то соответствующего значения случайной величины 2 нельзя найти в табличном задании закона распределения 2. В этом случае применяют следующую приближенную формулу:

,    (5.21)

основанную на том, что  оказывается асимптотически нормальным законом , z2 - есть 2% предел абсолютного уклонения нормальной переменной, заданный в табличных приложениях книг по теории вероятностей.

5.3. Формулы и алгоритмы для оценки результатов моделирования

При реализации моделирующего алгоритма на ЭВМ вырабатывается информация о состоянии моделируемых систем, которая представляет собой исходный материал для определения приближенных искомых величин по статистическим данным. Желательно так организовать фиксацию и алгоритмы обработки результатов моделирования, чтобы статистические оценки для искомых величин формировались постепенно по ходу моделирования, без специального запоминания всей информации о состояниях системы [8].

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

Р*(A)=m/N,       (5.22)

где m — число случаев (частота) наступления событий A, N - число реализаций (объем выборки). В данном случае для подсчета частоты достаточно предусмотреть один счетчик К, содержимое которого будет увеличиваться на единицу каждый раз при наступлении события A. Для получения значения Р*(A) после окончания моделирования содержимое счетчика К делится на N.

Если событие принимает значения в некоторой области величин, то область значений n случайной величины разбивается на отрезки так, что n={n1,n2…nm}, . Оценка вероятностей возможных i-х значений случайной величины определяется

Р*i(A)=mi/N,       (5.23)

где mi - число значений случайной величины в интервале ni. Для подсчета частоты необходимо предусмотреть m счетчиков К[I], содержимое которых будет увеличиваться на единицу каждый раз тогда, когда случайное событие A принимает значение из интервала ni.

Для получения значения Р*i(A) после окончания моделирования содержимое i-го счетчика К[I] делится на N. Алгоритм приведен на рис. 5.5. Примером непрерывной случайной величины A могут быть интервалы времени между движущимися автомобилями.

Если задать границы D(J), , где - заданное число границ оценки этой случайной величины А, то можно определить частоты событий А(J), состоящие в том, что значения случайной величины A меньше или равны границам D(J). Частоты А(J) записаны в счетчиках К(J), .

Рис. 5.5

Величина D(JМ) является наибольшей границей оценки случайной величины, т.е. D(1)<D(2)<…<D(JМ). Частота К(J) события А определена тем, что значение события меньше либо равно границе D(J).

На рис 5.6 приведен алгоритм подпрограммы STAT набора статистических данных. Входной переменной подпрограммы STAT является значение X непрерывной случайной величины A. В блоках 1, 2, 5 реализован цикл по переменной J. В блоке 3 проверяется условие, что значения X случайной величины A меньше или равны границам D(J). Если условие выполняется, то содержимое соответствующего счетчика К(J) увеличивается на единицу (см. блок 4).

В табл. 5.1 приведен пример частот некоторой случайной величины А. Обработка статистических данных, приведенных в табл. 5.1, позволяет построить кумулятивную эмпирическую функцию распределения.

Рис. 5.6

Таблица 5.1

Статистические данные результатов моделирования

Границы оценки

D(1)

2

D(2)

4

D(3)

6

D(4)

8

D(5)

10

D(6)

15

D(7)

20

D(8)

25

Номер счетчика

К(1)

К(2)

К(3)

К(4)

К(5)

К(6)

К(7)

К(8)

Частота события

37

100

193

240

280

362

425

500

Определяются частости появления события А в соответствии с формулами:

где Pj - теоретическое значение вероятностей.

Затем строится гистограмма кумулятивной эмпирической функции распределения по значениям Pj*. Пример построения приведен на рис. 5.7.

Рис. 5.7

Выдвигается гипотеза, состоящая в том, что найденная кумулятивная эмпирическая функция распределения может быть аппроксимирована известным теоретическим распределением P(x) (см. рис. 5.7). Проверка гипотезы осуществляется по критерию 2 (см. разд. 5.2.3).

Если определять частоты событий А(J), состоящие в том, что значения случайной величины A принадлежит интервалу (D(J+1)-D(J)), , и эти частоты записывать в счетчики К(J), , то алгоритм подпрограммы STAT в этом случае будет иметь вид, приведенный на рис. 5.8.

Рис. 5.8

Можно получить по статистическим данным  - частости попадания случайной величины A в интервалы (D(J+1)-D(J)), :

; ; ; …; .

Затем определить отношение  к величине j-го интервала (D(J+1)-D(J)):

Если при моделировании в счетчиках К(J) будут получены частоты событий, состоящих в том, что случайная величина А принадлежит интервалу (D(J+1)-D(J)), то частости  определятся:

Для построения кумулятивной эмпирической функции распределения частости  определятся следующим образом:

Среднее значение случайной величины определяется по формуле

,       (5.24)

где хk - возможные значения случайной величины, которые она принимает при различных реализациях процесса. На рис. 5.9 приведен алгоритм для определения среднего значения случайной величины.

Рис. 5.9

В этом алгоритме - такты моделирования; NZ - заданное число тактов моделирования; GEN(Х) – подпрограмма генерации случайной величины Х. После генерации всей выборки случайной величины Х в блоке 5 определяется среднее значение XS.

Оценкой S2* дисперсии случайной величины определится

,     (5.25)

где  - математическое ожидание случайной величины.

Эта формула неудобна, т.к. в процессе моделирования необходимо запоминать весь массив значений х1, х2, х3, …, хN. Известна упрощенная формула, согласно которой

,    (5.26)

т.е. для определения S2* достаточно в двух счетчиках накапливать значения  и . Для оценки корреляционного момента K случайных величин и с возможными значениями хk и yk применяется формула

.    (5.27)

Эта формула преобразуется к виду

,   (5.28)

требующему подсчета и запоминания в трех счетчиках соответствующих величин:

,    ,     .

Иногда искомыми величинами являются математическое ожидание и корреляционные функции случайного процесса Х(t). В теории случайных процессов изучаются закономерности изменения случайной величины от изменения неслучайного параметра, например времени, пространственной координаты и прочее. Основным понятием в теории вероятностей является понятие испытания с определенным множеством возможных элементарных событий - исходов испытания. Случайная величина X представляет однозначную числовую функцию X=f() элементарных событий, принимающего числовое значение в зависимости от исхода испытания.

Пусть каждому элементу множества соответствует не одно определенное значение, а определенная числовая функция f(t)(0,T) некоторого неслучайного параметра t. Так как для различных эти функции различны, то каждую такую функцию f(t) называют возможной реализацией случайного процесса Х(t). Совокупность всех возможных реализаций, т.е. множество функций f(t) образуют случайный процесс Х(t).

Распределение вероятностей случайного процесса Х(t) задают распределением вероятностей случайных величин Х(t1), Х(t2),, Х(ts), соответствующих любому конечному набору значений t1, t2,, ts параметра t (s=1,2,3,…).

На практике случайный процесс Х(t) определяют математическим ожиданием и дисперсией, являющимися функциями параметра t, а также корреляционной функцией. Рассмотрим, как определяют и как вычисляют эти функции. На рис. 5.10 показаны возможные реализации случайного процесса Х(t).

Математическим ожиданием случайного процесса Х(t) называется неслучайная функция МХ(t), значение которой при каждом значении t=ti равно математическому ожиданию МХ(ti) той случайной величины Х(ti), которая соответствует этому значению параметра.

Рис. 5.10

Математическое ожидание МХ(t) (см. рис. 5.10) представляет собой среднюю функцию, около которой группируются возможные реализации случайного процесса Х(t).

Дисперсией случайного процесса Х(t) называется неслучайная функция DХ(t), значение которой при каждом значении t=ti параметра t равно математическому ожиданию DХ(ti) той случайной величины Х(ti), которая соответствует значению параметра ti. Квадратный корень из дисперсии представляет среднее квадратичное отклонение случайного процесса Х(t) и определяется по формуле

.      (5.29)

Связь между случайными величинами Х(t*) и Х(t**), соответствующим значениям t* и t** случайного процесса Х(t), характеризуется их ковариацией

BX(t*,t**)=cov[Х(t*),Х(t**)]=

=M{[Х(t*)-MХ(t*)][Х(t**)-MХ(t**)]}.   (5.30)

Ковариация представляет собой неслучайную функцию BX(t*,t**) двух переменных t* и t**, которая графически может быть представлена поверхностью, как это показано на рис. 5.11.

Рис. 5.11

Функция BX(t*,t**) называется корреляционной функцией или автокорреляционной функцией случайного процесса Х(t).

Интересующий интервал (0,T) разбивается на части с шагом t. Накапливают значения Хk(ti) реализаций случайного процесса Х(t) для фиксированных моментов времени ti. Затем вычисляют оценки для математического ожидания по формуле

.     (5.29)

Оценки для корреляционной функции BX(t*,t**) вычисляются по формуле

,   (5.30)

где t* и t** «пробегают» все значения t. Так при моделировании при применении формулы (5.30) необходимо накапливать N значений Xk(t*) и N значений Xk(t**). На практике для оценки корреляционной функции BX(t*,t**) применяют формулу

.    (5.31)

При применении формулы (5.31) необходимо три счетчика для подсчета сумм

, , .


6. МОДЕЛИРОВАНИЕ ВЕРОЯТНОСТНЫХ АВТОМАТОВ

6.1. Аналитическое определение вероятностных автоматов

6.1.1. Формальное задание и классификация. Вероятностные автоматы (ВА) относятся к дискретно-стохастическому классу моделей. Данный тип моделей служит инструментом изучения динамических систем, имеющих стохастическую природу функционирования с дискретным временем. ВА является типичным представителем таких систем (probabilistic automat) и носит название P-схемы или P-автомата. В общем виде такой автомат является потактным преобразователем информации с памятью, функционирование которого может быть описано статистически.

Математический аппарат ВА применим для разработки методов проектирования дискретных систем, проявляющих статистически закономерное случайное поведение, для выяснения возможностей таких систем и обоснования границ целесообразности их использования, а также при решении различных задач синтеза. Аппарат ВА применяется также для моделирования дискретно-стохастических объектов, у которых подача входных параметров, изменение состояния и формирование выходных параметров осуществляется в дискретные моменты времени ti (t0,t1,...,ti…). Состояние объекта определяется через предшествующие состояния и входной параметр. Выходной параметр определяется через состояние в данном такте времени, состояние в предшествующем такте, а также через входной параметр.

Для формального описания ВА следует задать распределение начальных состояний, множество входных параметров Х={х12,...,хm}, множество состояний Z={z1,z2,...,zn}, множество выходных параметров Y={y1,y2,...,yr}. Элементы множества Х,Z,Y называют входным, внутренним и выходным алфавитом.

Определение. Вероятностным автоматом называется математическая схема, которая задается следующим набором [7,14]:

ВА=<Z,Y,Р0,{Р(zt,yt/zt-1t)}>,

где Р0 - распределение начальных состояний, Р0=||||, - вероятность того, что в такте времени t0 автомат будет находиться в состоянии zi; Р=||Р(zt,yt/zt-1t)|| - стохастическая матрица, в которой Р(zt,yt/zt-1t)=Р{z(t)=zt, y(t)=yt/z(t-1)=zt-1, х(t)=хt} - условная вероятность того, что в такте времени t автомат будет в состоянии zt, на выходе будет иметь параметр yt при условии, что в такте t-1 автомат был в состоянии zt-1, а на вход был подан параметр хt. При моделировании следует определить функции переходов и выходов. Функцию переходов задают в виде стохастической матрицы ||Р{zt(t)=z(t)/zt-1t}||.

Функция выходов определяет выходные параметры и задается в виде стохастической матрицы ||Р(yt/zt-1t,zt)||, в которой Р(yt/zt-1t,zt)=Р{y(t)=yt/z(t-1)=zt-1,х(t)=хt, z(t)=zt}.

Определим условную вероятность Р(yt/zt-1tzt):

Р(zt,yt/zt-1t)=Р(zt/zt-1t)Р(yt/zt-1t,zt).

Просуммируем правую и левую части по всем значениям yi и получим

Сумма в правой части равна единице, так как это сумма вероятностей полной группы событий. Тогда вероятность Р(yt/zt-1t,zt) определится формулой

.

6.1.2. Классификация ВА. Классификация ВА зависит от способов определения вероятности Р(yt/zt-1t,zt) функции выходов и вероятности Р(yt/zt-1t) функции переходов.

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

Р(yt/zt-1t,zt)=Р(yt/zt-1t), (автомат Мили).

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

Р(yt/zt-1t,zt)=Р(ytt,zt).

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

Р(yt/zt-1t,zt)=Р(yt/zt-1,zt).

Существует правильный ВА первого рода, у которого

Р(yt/zt-1t,zt)=Р(yt/zt-1),

и правильный вероятностный автомат второго рода, у которого

Р(yt,zt-1t,zt)=Р(yt/zt), (автомат Мура).

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

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

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

Вероятностный автомат первого рода с детерминированной функцией выходов называется марковским.

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

ВА=<Х,Z,Р0{Р(zt/zt-1t}>.

Если мощность множества Z равна единице, то такой автомат называется автоматом без памяти.

Если мощность множества Х равна единице, то такой автомат называется автономным.

Автономный абстрактный ВА называется дискретной цепью Маркова и задается в следующем виде:

ВА=<Z,Р0{Р(zt/zt-1)}>.

6.2. Табличное задание функций переходов и выходов

Задание условных вероятностных мер Р(zt,yt/zt-1t) возможно как задание стохастического отображения ZХZY табличным способом. В табл. 6.1 приведен вид совместного задания функций переходов и выходов.

Таблица 6.1

Совместное задание функций переходов и выходов

ZХ

ZY

z1y1

z1y2

z1yr

zny1

zny2

znyr

z1х1

z1хm