Образовательный стандарт дисциплины "Системное моделирование"

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

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

Возникновение таких ПЭВМ позволило колоссально расширить интеллектуальные возможности человека, раскрепостить его умственные силы.

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

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

Вопрос 1

Что понимается под названием вычислительный эксперимент?

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

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

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

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

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

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

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

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

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

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

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

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

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

Таким образом, можно выделить три основные составные части метода статистических испытаний:

построение математической модели и моделирующего алгоритма исследуемой системы;

2) формирование случайных величин с заданным законом распределения вероятностей;

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

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

Статистическая обработка и оценка точности результатов моделирования основывается на предельных теоремах теории вероятностей: теореме Чебышева и теореме Бернулли.

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

, (1)

где - сколь угодно малое положительное число,

.

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

. (2)

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

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

. (3)

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

.

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

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

, (4)

где - функция Лапласа (интеграл вероятностей);

- аргумент функции Лапласа;

- среднее квадратическое отклонение величины .

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

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

. (5)

Положим

, (6)

тогда получим

. (7)

Сравнивая выражения (3) и (7), найдем условие, при котором ошибка метода не превысит величину с вероятностью :

. (8)

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

. (9)

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

. (10)

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

Задавая доверительную вероятность , получаем из формул (5), (8) доверительную оценку

(11)

с надежностью . Отсюда вытекает, что ошибка метода статистических испытаний пропорциональна величине . Следовательно, чтобы уменьшить ошибку в 10 раз (то есть, чтобы получить в ответе еще один верный знак), нужно увеличить число испытаний в 100 раз. Чтобы получить достаточно высокую точность, необходимо провести тысячи испытаний. Метод особенно эффективен при решении задач, в которых результат нужен с точностью порядка 5 - 10 %.

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

Так как есть искомая величина, а - ее приближенное значение, то есть ошибка метода.

Введя снова величину , удовлетворяющую неравенство , и доверительную вероятность , получим на основании теоремы Бернулли

. (12)

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

, (13)

где также находится из условия (8).

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

Основными достоинствами метода статистических испытаний являются:

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

простота осуществимости. Составляется программа для одного испытания, затем испытание повторяется раз. Нет необходимости в создании специальных устройств;

простота оценки точности полученных результатов.

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

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

2. Моделирование системы массового обслуживания

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

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

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

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

,

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

- случайная величина, равномерно распределенная в интервале .

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

.

Если в момент поступления заявки все каналы заняты, система выдает отказ.

-й канал

Рис.2.1 Показатели эффективности системы за время функционирования

Требуется определить показатели эффективности системы за время функционирования .

Рис. 2.2 Логическая схема алгоритма процесса обслуживания заявок

Введем следующие обозначения (рис. 2.1):

- момент поступления -й заявки;

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

- момент освобождения -го канала:

За начальный примем момент поступления первой заявки .

В этот момент все каналы свободны .

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

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

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

Оператор 3 определяет, принадлежит ли -я заявка заданному интервалу времени функционирования системы. Если условие выполняется, заявка поступает на обслуживание, управление передается оператору 4. В противном случае испытание заканчивается и управление передается на счетчик числа испытаний.

Оператор 4 сравнивает между собой моменты освобождения каналов системы и выбирает канал, освободившийся раньше всех. Если есть каналы, освободившиеся одновременно, выбирается канал с меньшим номером. Пусть номер выбранного канала равен . Заявка поступает для обслуживания в этот канал. Очевидно, в каждом испытании первая заявка поступает в первый канал, вторая заявка - во второй канал, третья - в третий, так как в начальный момент .

Оператор 5 сравнивает момент освобождения выбранного канала с моментом поступления -й заявки . Если условие не выполняется, это означает, что все каналы к моменту заняты, система выдает отказ и управление передается оператору 6. Если же указанное условие выполняется, канал с номером свободен и управление передается оператору 10.

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

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

Оператор 7 формирует поток заявок, то есть по формуле

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

Оператор 8 формирует момент поступления следующей заявки

.

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

Если условие , проверку которого осуществляет оператор 5, выполняется, то выбранный оператором 4 канал с номером свободен и заявка передается в этот канал для обслуживания. Дальше необходимо определить время обслуживания -й заявки, найти время освобождения канала, зафиксировать обслуженную заявку и формировать следующую заявку. Эту задачу выполняют операторы 10, 11 и 12.

Оператор 10 определяет время обслуживания -й заявки согласно формуле

.

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

.

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

С оператора 12 управление передается на оператор 7 и дальше формируется следующая заявка так же, как и в рассмотрением случае отказа в обслуживании.

Если неравенство не выполняется (следовательно , это означает, что -я заявка уже не принадлежит заданному интервалу, и реализация на этом заканчивается.

Оператор 13 представляет собой счетчик числа испытаний.

Оператор 14 проверяет, получено ли уже заданное число испытаний . Если неравенство выполняется, управление передается оператору 15.

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

Если неравенство не выполняется, управление передается оператору 16.

Оператор 16 осуществляет статистическую обработку полученных результатов и вычисляет требуемые показатели эффективности функционирования системы за время .

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

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

Пусть мы хотим моделировать систему массового обслуживания, поступление требований в которой подчинено пуассоновскому распределению со средним 3 клиента в час, а время обслуживания равно 0,2 ч с вероятностью 0,5 или 0,6 ч с вероятностью 0,5. Клиенты обслуживаются согласно дисциплине «первым пришел - первым обслуживаешься»; длина очереди, а также источник поступления клиентов не ограничены. Предположим, что в начальный момент моделирования клиентов нет.

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

.

Поскольку время обслуживания равно либо 0,2, либо 0,6 ч с равными вероятностями, время обслуживания определяется как

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

Событие, связанное с поступлением клиента

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

Проверка состояния системы (простой или работа).

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

б) Если система работает, поставить поступившего клиента в очередь и увеличить ее длину на единицу.

Событие, связанное с окончанием обслуживания

Проверка состояния очереди (пустая или непустая).

а) Если очередь пуста, объявить простой системы.

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

ч.

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

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



Поступление Поступление Конец обслуживания

Рис. 2.3 События, связанные с окончанием обслуживания

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

.

Система объявляется работающей, а время простоя корректируется следующим образом: Время простоя ч.

Осуществившиеся до настоящего момента события показаны на рис. 2.3.

Следующее по времени событие - поступление требования в момент . Поскольку система продолжает работать, требование ставится в очередь, а длина очереди корректируется:

(в момент ).

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

.

(В рассматриваемом примере полезно наносить новые события на рис. 2.3 по мере их получения.)

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

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

(в момент ),

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

.

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

(в момент ),

а суммарное время ожидания становится равным

ч.

Доля времени простоя системы, %	=	Суммарное время простоя Период моделирования	100
Среднее время ожидания клиентом обслуживания	=	Суммарное время ожидания Число поступающих клиентов

Используя , получаем время завершения обслуживания данного клиента:

.

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

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

Средняя длина очереди =	Площадь Моделируемый период

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

Средняя длина очереди = A/t

Площадь A

1,07 1,09 1,14 Моделируемое время t

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

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

Процедура имеет большое сходство с физическим экспериментом.

Выполнив упражнение с помощью «ручных» вычислений, можно убедиться в необходимости использования ПЭВМ при моделировании. Использование ПЭВМ становится еще более привлекательным из-за возможности применения таких специализированных языков моделирования, как GASP, SIAM, GPSS и SIMSCRIPT. Эти языки разработаны для того, чтобы избавить пользователя от утомительной необходимости программирования многочисленных деталей. Например, все языки дают возможность автоматически генерировать и запоминать события в хронологическом порядке с помощью всего одного оператора. Кроме того, все языки обладают очень простыми операторами для автоматического табулирования операционных характеристик системы. Имея подобные языки, пользователь может сосредоточить усилия на улучшении модели.

3. Получение результатов наблюдений при моделировании

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

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

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

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

1) очень длительных прогонов модели, позволяющих увеличить вероятность достижения стационарного состояния;

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

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

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

затраты на моделирование могут существенно зависеть от продолжительности прогонов модели;

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

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

Рассмотрим два метода получения наблюдений: 1)метод повторения, 2)метод подынтервалов. Имеются и другие методы, однако эти два, по-видимому, наиболее подходят для практических приложений.

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

Метод повторения

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

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

и -й доверительный интервал для точного значения

среднего вычисляется как

,

и есть -статистика с степенями свободы.

Метод подынтервалов

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

Рис. 3.1 Графическая интерпретация метода подынтервалов

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

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

4. Прикладные задачи имитационного моделирования

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

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

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

Значение изменяющегося параметра объекта прогнозирования для каждого момента на периоде основания можно представить в виде

,

где - значение динамического ряда в -й момент времени (год) периода основания;

- значение динамического ряда в предыдущий момент времени;

- приращение переменной объекта прогнозирования в -й момент времени по сравнению с предыдущими;

- число значений динамического ряда.

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

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

,

где - номер шага на периоде упреждения;

- число шагов на периоде упреждения;

- значение переменной объекта прогнозирования на предыдущем шаге;

- моделируемое значение приращения на -м шаге.

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

,

где - точечный прогноз на -й период упреждения;

- конечное значение динамического ряда.

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

; (4.1)

, (4.2)

где - число реализаций точечного прогноза;

- разыгрываемое значение приращения на -м шаге периода упреждения в -й реализации точечного прогноза;

- значение -й реализации точечного прогноза, определяемое по зависимости (1).



Рис. 4.1 Графическое отображение процесса случайного блуждания

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

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

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

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

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

,

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

Рис. 4.2 График определения предыдущих и последующих приращений

Графически определение предыдущих и последующих приращений показано на рис. 4.2.

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

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

.

Множество последующих приращений определяется по зависимости

или

.

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

(4.3)

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

 (4.4)

Оценка значения коэффициента корреляции определится по зависимости

. (4.5)

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

(4.6)

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

При моделировании приращений на следующих шагах периода упреждения

.

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

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

; (7)

,

где - значение выборочного коэффициента корреляции, определяемое по зависимости (4.5).

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

, (4.8)

где - нормированная нормально распределенная случайная величина, моделируемая с помощью алгоритма.

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

. (4.9)

Рис. 4.3. Блок-схема алгоритма прогнозирования с использованием ориентированного процесса случайного блуждания

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

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

вычисление случайного значения по зависимостям (4.8) и (4.9);

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

вычисление приращения по зависимости (4.6) при полученном в п. 2 значении коэффициента корреляции , определенном в п. 3 значении .

Многократно имитируя приращения и используя зависимости (4.1) и (4.2), вычисляются характеристики прогноза. Блок-схема алгоритма изображена на рис.4.3.

К достоинствам рассмотренного метода прогнозирования относятся:

простота вычислительного алгоритма;

возможность использования при ограниченной на периоде основания информации (начиная с 7-9 значений динамического ряда);

простота оценивания точности прогноза (определения дисперсии).

4.2 Модифицированный имитационным моделированием метод экспоненциального сглаживания

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

, (4.10)

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

,(4.11)

аппроксимирующего выражение (4.10), несколькими членами .

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

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

, (4.12)

где ;

- погрешность, обусловленная приближенностью исходной информации;

- погрешность, связанная с методом прогнозирования;

- погрешность, вызванная неточностью вычислений;

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

Одной из наиболее весомых является методическая ошибка, зависящая от числа членов разложения. В работах [1], [2] приводятся аналитические зависимости для выполнения параметров аппроксимирующего многочлена при . Вывод таких зависимостей для представляет значительные трудности. Кроме того, любое увеличение числа членов выражения (4.11) влечет за собой потребность увеличения объема исходных данных, необходимых для определения оценок начальных значений коэффициентов (методом наименьших квадратов или в более общем случае методом максимального правдоподобия), далее предлагается модификация метода экспоненциального сглаживания, основанная на принципах группового учета аргументов. Сущность метода заключается в том, что математическая модель объекта прогнозирования

,

называемая в соответствии с терминологией работы [1] его «полным описанием», заменяется набором «частных описаний» вида

.

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

.

Аналогичная процедура повторяется до тех пор, пока величина критерия фильтрации уменьшается или увеличивается в зависимости от его содержания (при этом исходная информация делится на две выборки: обучающую и проверочную). Для практических расчетов в качестве такого критерия рекомендуется принимать среднеквадратическую ошибку аппроксимации модели на проверочной выборке, которая, как установлено в работе [10], при увеличении числа уровней фильтрации, а, следовательно, сложности модели, достигает экстремального значения. Сложность модели (измеряется числом ее членов), соответствующая экстремальному значению критерия, является оптимальной. На последнем уровне фильтрации фиксируется «частное описание», значение которого минимально. На предпоследнем уровне выбираются «частные описания», являющиеся аргументами последнего уровня, и т.д. Так как «частные описания» являются функцией двух аргументов, их коэффициенты легко определяются по небольшому количеству исходных данных. Исключая промежуточные переменные можно получить модель исследуемых характеристик объекта прогнозирования в виде аналога «полного описания»

,

где в общем случае .

Как известно, особые трудности при увеличении числа членов в разложении Тейлора связаны с получением аналитических зависимостей для определения вектора коэффициентов . Из работы [2] следует, что

,

где - вектор-столбец размером сглаженных значений процесса

;

- вектор-столбец размером неизвестных коэффициентов

;

- матрица размером , элементы которой, соответствующие -й строке и -му столбцу, вычисляются по зависимости

. (4.13)

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

вектор выражается зависимостью . (4.14)

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

, (4.15)

где ;

Рис. 4.4 Блок-схема алгоритма прогнозирования по методу модифицированного экспоненциального сглаживания

- коэффициенты многочлена с переменной .

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

, (4.16)

где .

Рис. 4.4. Блок-схема алгоритма прогнозирования по методу модифицированного экспоненциального сглаживания (продолжение)

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

.

Блок-схема алгоритма прогнозирования, составленного в соответствии с изложенными положениями, изображена на рис.4.4.

Автоматический подбор вида экстраполируемой функции

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

Известно [1], [2], [3], что любой процесс можно представить в

, (4.17)

где - исходный процесс (функция одного переменного);

- приближенная модель процесса (описание с помощью ФГС);

- остаток (некоторая функция точности приближения).

В наиболее общем виде ФГС для одного аргумента записывается в виде [1], [2]

, (4.18)

где - некоторое фиксированное натуральное число;

- начальное значение фактора-аргумента на рассматриваемом интервале;

- постоянные действительные параметры;

- специальный (степенной) определитель -го порядка;

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

, .

При функция с гибкой структурой имеет вид

, (4.19)

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