CQ
О
о
о
о
с*
S
и
ш
Я. Я, Бусленко
МЕТОД
СТАТИСТИЧЕСКОГО
МОДЕЛИРОВАНИЯ
СТАТИСТИКА МОСКВА 1970

Брошюра посвящена методу
статистического моделирования,
реализуемому на быстродействующих
электронных цифровых машинах.
Основное внимание обращено на методы
построения машинных моделей для
различных сложных систем,
встречающихся в народном хозяйстве,
организации производства, учете,
сборе и обработке экономической
информации.
Динамика функционирования
сложных систем иллюстрируется
примерами систем массового
обслуживания различных типов и близких к
ним формализованных схем. В связи
с этим в брошюре излагаются
методы формирования потоков
однородных событий, имитации работы
обслуживающих приборов с учетом
множества действующих факторов, а
также методы построения
алгоритмов для моделирования сложных
систем на ЭВМ.
Брошюра предназначена для
инженеров-экономистов, экономистов и
статистиков, интересующихся методами
машинного моделирования, а также
для студентов и аспирантов
экономических вузов, связанных с
применением средств вычислительной техники.
Редколлегия серии
«Математическая статистика
для экономистов»
Боярский А. #., Венвцкий И. Г., Го-
ленко Д. #., Дружинин Н. /С, Пасха-
вер И. С., Привезенцева А. Г.
1—8—3
22—70

ВВЕДЕНИЕ
Социалистическое народное хозяйство
развивается высокими темпами. С каждым годом усложняется
управление экономикой нашей страны. Оно
превращается в важнейший объект механизации и автоматизации,
становится исключительно широкой сферой для
реализации новейших научных идей, для применения
вычислительной техники и различных
информационно-технических комплексов.
Очень остро подчеркнул значение этой проблемы
декабрьский A969 г.) Пленум ЦК КПСС. На Пленуме
отмечалось, что современные условия требуют
совершенствования методов, самой системы управления, в
частности работы, связанной со сбором, быстрой обработкой
и анализом информации, дальнейшей разработки
научно обоснованных методов принятия решений,
существенных поправок в структуре управления, что позволяет
строить ее более рационально, избавляясь от ненужных
звеньев, устанавливая правильное соотношение прав и
обязанностей, власти и ответственности на всех
уровнях1.
Основу применения электронной вычислительной
техники в области переработки плановой и учетной
информации, в системах управления экономическими и
производственными процессами составляют
математические методы. К настоящему времени широкие круги
экономистов, работников учета, статистики уже хорошо
познакомились с методологическими особенностями
решения экономических задач при помощи математических
методов. Широко освещены в литературе идеи и
практическая реализация линейного программирования,
построения сетевых графиков. Однако этим далеко не ис-
См. «Правда» от 13 января J970 г,

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

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

Большое значение приобретает метод
статистического моделирования для обработки данных, получаемых
путем обследования состояний сложной системы и
экспериментального изучения процессов ее
функционирования. Существо этого вопроса состоит в следующем. В
последнее время обследованиям подвергаются весьма
сложные системы, функционирование которых зависит
от большого числа параметров. К сожалению, во многих
случаях некоторые параметры системы могут быть
оценены лишь косвенно, по результатам их влияния на
характеристики системы, непосредственная регистрация
которых не представляет затруднений. Очевидно, что для
оценки такого рода параметров необходимы
соотношения между регистрируемыми и оцениваемыми
величинами. Эти соотношения могут быть получены в результате
статистического моделирования системы.
Приемы машинной обработки данных с
использованием статистических моделей в настоящее время
интенсивно разрабатываются.
Глава V посвящена сравнительно новой проблеме —
моделированию автомобильного движения — и является
в известном смысле • практической иллюстрацией
применения метода статистического моделирования для
решения весьма сложных задач. Она написана В. Н. Бус-
ленко по материалам, разработанным совместно с
Э. В. Кручининой.

Глава I
МЕТОД СТАТИСТИЧЕСКИХ ИСПЫТАНИЙ
(МОНТЕ КАРЛО)
1. Сущность метода
Одним из численных методов, получивших
распространение при появлении быстродействующих
электронных вычислительных машин (ЭВМ), является метод
статистических испытаний, или метод Монте Карло. Он
базируется на использовании так называемых случайных
чисел — возможных значений 'некоторой случайной
величины с заданным распределением вероятностей.
При реализации метода статистических испытаний на
ЭВМ случайные числа вырабатываются специальной
электронной приставкой к ЭВМ (датчиком случайных
чисел) или самой машиной по специальной программе.
В любом из этих вариантов для генерирования
случайных чисел используется аппаратурная (датчик
случайных чисел) или алгоритмическая (специальная
программа для ЭВМ) модель некоторого случайного процесса,
вероятностные характеристики которого известны или
могут быть оценены экспериментально. Ранее, до
появления ЭВМ, этой цели служили простейшие случайные
процессы, такие, как бросание монеты (выпадение
герба или решки с вероятностями pi ==р2 = 0,5), бросание
игральной кости (/?;= -g~; /=1,2,..., 6), извлечение карт
из тщательно перетасованной .колоды (Рг= ТШ"• *=1>2'
..., 36), вращение рулетки и т. д.1.
1 Между прочим, отсюда происходит и наименование «метод
Монте Карло» — по названию столицы княжества Монако на берегу
Средиземного моря, известного игорными /домами.

Сущность метода статистических испытаний поясним
на примерах.
В качестве первого примера рассмотрим вычисление
площади некоторой фигуры произвольной формы.
Остановимся сначала на частном случае решения этой
задачи. Пусть требуется вычислить площадь фигуры (см.
рис. 1), ограниченной отрезками О А и 01 на осях пря-
Т
Рис. I.
моугольных координат ОХ и 0Y, кривой y=f (х) и
ординатой IS (искомая площадь на рис. 1 заштрихована),
причем будем считать выполненным условие O^f(x)^l
для всех х\ O^jc^I.
Пользуясь обычными численными методами для
приближенного вычисления искомой площади, поступают
следующим образом. Разбивают отрезок @,1) на оси
ОХ на п равных частей длиной Ддг= —. Искомую
площадь 5 представляют в виде суммы площадей п
элементарных фигур, ка.к показано на рис. 2. Площадь каждой

элементарной фигуры приближенно можно заменить
площадью соответствующего прямоугольника, равной AS* =
= А* •/(**), где Хг — некоторая точка на оси ОХ внутри
/-го интервала. Точки х\ выбирают таким образом, чтобы
площадь -прямоугольника была возможно более близкой
к площади элементарной фигуры. Очевидно, что при до-
У
О *
Рис. 2.
статочно большом п точность вычисления площади
можно сделать вполне приемлемой. Отметим для
дальнейшего, что такой способ определения площади требует
вычисления значений функции f(x) в п точках.
Посмотрим теперь, как решается эта же задача
методом статистических испытаний. Пусть мы имеем
случайную величину ?, равномерно распределенную на
отрезке [0, 1] (см. [5]). Это значит, что вероятность
попадания ее возможных значений х в интервал (а, Ь)
пропорциональна длине интервала и не зависит от
местоположения его на отрезке [0, 1].

Если возможные значения х равномерно
распределенной случайной величины | заполняют отрезок [0, 1]
на оси ОХ и возможные значения у случайной величины
г] заполняют тот же отрезок на оси OY, то пары чисел
(х, у) определяют случайную точку (х, у) на плоскости
OXY, имеющую равномерное распределение в квадрате
@,0), @,1), A,1), A,0), который мы в дальнейшем
будем называть единичным квадратом. Это значит, что
вероятность попадания точки (х,у) в некоторую область
а пропорциональна площади этой области и не зависит
от расположения ее внутри единичного квадрата.
Рис. 3.
Проведем мысленный эксперимент: внутрь
единичного квадрата случайным образом с равномерным
распределением бросается точка. Это эквивалентно выборке
пары чисел xt и уи являющихся возможными
значениями ? и г] соответственно. После N таких испытаний (где
N достаточно велико) на плоскости появится N
случайно расположенных точек, равномерно распределенных в
единичном квадрате (см. рис. 3).
10

Предположим, что количество точек под кривой
y = f(x) равно га, а над кривой y = f(x) равно N—m
(точки, попадающие точно на кривую, будем считать
находящимися под кривой).
Если следовать геометрическим соображениям, ясно,
что вероятность Р попадания точки в часть квадрата,
находящуюся под кривой # = / (х), равна отношению
площади S этой части квадрата к площади всего
квадрата. Частота -^ попадания точки в часть квадрата под
кривой y — f (x) при достаточно большом N близка к
вероятности Р (см. [5]). Отсюда следует, что в качестве
приближенного значения искомой площади можно взять
т
частоту -др-, т. е.
т
S» —• A.1)
N '
Для решения рассмотренного примера на ЭВМ «нет
необходимости в воспроизведении всех указанных выше
действий. Сущность метода статистических испытаний
для данного случая состоит в моделировании
эксперимента при помощи случайных чисел.
Процедура решения выглядит следующим образом:
1. Выбирается случайное число & из отрезка [0; 1] с
равномерным законом распределения (из таблиц
случайных чисел или вырабатывается самой машиной с
помощью датчика случайных чисел); это случайное число
принимается в качестве координаты случайной точки Хг
на оси ОХ.
2. Вычисляется значение рассматриваемой функции
f (Xj) в точке Xj.
3. Вырабатывается следующее случайное число |г+ь
принимаемое в качестве координаты точки у$ на оси О У;
таким образом Xj=^i и #j = Im-i определяют случайную
точку на плоскости внутри единичного квадрата.
4. Количество выработанных таким образом
случайных точек (пар случайных чисел) подсчитывается
специальным счетчиком, который мы будем называть
счетчиком количества испытаний N.
5. Значение функции / (х^ сравнивается со
случайным числом |,-+1. Если неравенство (см. рис. 3)
E<+i<f(*j) A-2)
и

выполнено, что соответствует попаданию случайной
точки (Xj, t/j) в часть квадрата под кривой y = f{x), то
результату сравнения присваивается специальный признак
@=1, если не выполнено (о = 0.
6. Полученные значения признака со прибавляются
к содержимому счетчика количества точек под кривой
(т).
7. Управление передается снова первой операции,что
соответствует переходу к новой случайной точке (*j+i;
i/j+i). После проведения N таких экспериментов
определяется приближенное значение площади под кривой:
S=--
N
Рассмотренная процедура не требует запоминания
всех случайных чисел, полученных в результате
эксперимента. Запоминаются только значения т и N. Это
немаловажное обстоятельство вообще характерно для
реализации метода статистических испытаний на ЭВМ.
Точность решения задачи методом статистических
испытаний растет с увеличением количества испытаний
N и при достаточно больших N становится приемлемой
с практической точки зрения. (Этот вопрос будет более
обстоятельно выяснен ниже, в § 5.)
Целесообразно обратить внимание на ряд возможных
обобщений.
Во-первых, метод статистических испытаний
позволяет вычислять площади фигур произвольной формы и
любых размеров. Пусть, например, требуется вычислить
площадь S' фигуры, изображенной на рис. 4. Поскольку
прямоугольник АВЕД не является единичным,
целесообразно изменить масштаб по осям координат:
Х г У
ДВ * АЕ
Тогда квадрат АВЕД будет единичным, а искомая
площадь
-AE-S*,
где 5* — площадь фигуры, выраженная в единицах
измерения, соответствующих новым масштабам.
12

Величина S* может быть определена методом
статистических испытаний. Процедура решения задачи в
основном совпадает с рассмотренной выше, за
исключением того, что теперь для каждого xj нужно вычислять два
значения у'{х}) и у(х,) ординат точек, лежащих на
границе фигуры (см. рис. 4), и проверять
справедливость неравенства
:х<). A.3)
причем уг выбирается уже не из отрезка [0,1], а из
преобразованного (сдвинутого) единичного отрезка.
Если это неравенство выполнено, случайная точка
(х-, иЛ находится внутри контура фигуры и тогда ш=1.
Во-вторых, задача вычисления площади является
частным случаем более общей задачи интегрального
исчисления. В самом деле, площадь 6 (см. рис. 6)
может быть выражена как
13

i
S f(x)dx9 A.4)
поэтому процедура определения площади
одновременно является процедурой вычисления интеграла вида
A.4). Напомним, что здесь
Если же это условие не выполнено и пределы
интегрирования произвольны, необходимо преобразовать
масштабы по осям координат. Например, пусть
требуется вычислить интеграл
ь
/= J f(x)dxt A.5)
а
где максимальное значение f(x) в отрезке [а,Ь] равно
fx. Используя замену переменных
x = a+(b-a)z A.6)
и изменение масштаба по оси у, получим
ь 1
/= $f(x)dx = fmax(b-a) §f*(z)dz, A.7)
а О
где
^a)*]. A.8)
Таким образом, рассматриваемая задача сводитсяс к
предыдущей, решаемой с помощью интеграла вида A.4).
2. Попадание случайной величины в заданный интервал.
Среднее значение функции /от случайной величины
В дальнейшем мы будем различать дискретные и
непрерывные случайные величины. Дискретная
случайная величина | может принимать лишь конечное или
счетное множество возможных значений Х{. Каждому
14

возможному значению Xi ставится в соответствие его
вероятность р*. Зависимость вида
будем называть законом распределения дискретной
случайной величины. При исследовании дискретных
случайных величин пользуются также функцией
распределения F(x), которая при каждом х равна вероятности
того, что |<*
F(x)=P(l<x).
Если известен закон распределения Pi^Pi{Xi),
всегда можно определить функцию распределения
х{<х
и, наоборот, 'по заданной функции распределения F(x)
можно найти закон распределения
Кроме дискретных случайных величин
рассматриваются случайные величины, возможные значения которых
сплошь заполняют соответствующие интервалы на
числовой оси. Такая случайная величина характеризуется
функцией распределения, которая определяется также,
как и для дискретных случайных величин; при каждом
х функция распределения
F(x)=P(l<x).
Во многих случаях функция распределения F(x)
оказывается дифференцируемой. Если существует
производная
случайная величина | называется непрерывной, а
функция f(x)—функцией плотности вероятностей случайной
величины ?. Очевидно, что
X
= lf(x)dx.
— со
15

Рассмотрим случайную величину |, возможные
значения которой х принадлежат некоторому интервалу на
оси ОХ. Закон распределения этой случайной величины
задан функцией плотности f(x). Вычислим вероятность
попадания случайной величины внутрь полуинтервала
[а, Ь), целиком лежащего в вышеуказанном интервале
на оси ОХ.
Из теории вероятностей известно (см. [5]), что
ь
f(x)dx.
Для решения этой задачи воспользуемся методом
статистических испытаний.
Легко видеть, что в данном случае применима уже
изложенная нами процедура вычисления интеграла вида
A.9). Однако мы рассмотрим другой вариант метода
статистических испытаний. Предположим, что в нашем
распоряжении имеются случайные числа Хь закон
распределения которых описывается функцией плотности
/ (х). Поскольку случайные числа х\ можно
рассматривать как возможные значения случайной величины |,
искомая вероятность P(a^g<b) A.9)) приближенно
(при достаточно большом числе испытаний N) равна
частоте попадания случайных чисел х\ в полуинтервал
а^хг<Ь.
Р&-, A.10)
N
где т — количество случайных чисел, удовлетворяющих
неравенству
а^КЬ A.11)
при ./V испытаниях.
Другими словами, в условиях данной задачи
ъ
IWdxKj- A.12)
а
Процедура вычисления интеграла на ЭВМ
следующая:
1. Из множества случайных чисел с законом
распределения f (x) выбираем значение х*.
16

2. Xi сравниваем с граничными значениями
полуинтервала (а,Ь); выполнение неравенства a^Xi<b
отмечаем признаком (о=1, невыполнение — о) = 0.
3. Полученное значение признака со прибавляем к
количеству случайных чисел, попавших внутрь
полуинтервала (а, Ь).
4. К содержимому счетчика количества испытаний
N прибавляется единица.
5. Управление передается снова первой операции.
Проведя /V таких испытаний, вычисляем значение
частоты -jjr.
Нетрудно заметить, что рассмотренный второй
вариант подхода к вычислению интеграла
существенным образом отличается от первого варианта
(см. стр. 14) подхода. Так, в первом варианте для
реализации каждого испытания были необходимы два
случайных числа Х{ и уь Кроме того, требовалось
вычисление подынтегральной функции для каждого значения х\
случайной величины |, что, вообще говоря, может быть
связано со значительным объемом вычислений. С
другой стороны, первый способ является более общим
случаем вычисления интеграла; этот способ пригоден для
подынтегральных функций весьма широкого класса,
тогда как второй способ предполагает, что
подынтегральная функция является функцией плотности некоторой
случайной величины, что накладывает на нее
определенные ограничения. Первый способ оперирует с равномерно
распределенной случайной величиной, которую легко
получить на ЭВМ, тогда 'как рассмотренный здесь требует
специального 'преобразования равномерного
распределения в заданное.
Перечисленные отличия разграничивают классы
задач, которые удобно решать тем или другим способом.
Так, первым способом широко пользуются для
вычисления многократных интегралов, в то время как второй
способ находит применение при моделировании
сложных систем методом статистических испытаний.
17

Перейдем к рассмотрению еще одной разновидности
применения метода статистических испытаний для
вычисления интегралов.
Пусть ?— случайная величина, возможные значения
которой принадлежат отрезку [a/b], a h(x)—функция
плотности этой случайной величины. Рассмотрим
непрерывную функцию r\ = g(?)- Из теории вероятностей
известно (см. [5]), что среднее значение, или
математическое ожидание функции g (|), равно
M[g(l))= $ g(x)k(x)dx. A.13)
а
Воспользуемся методом статистических испытаний
для вычисления интегралов вида A.13). Из
совокупности возможных значений х случайной величины §
выберем последовательность Х\\ х2\ ...; Xn и вычислим
значения g(Xi); i'=l, 2, ..., N; причем хи не
принадлежащие отрезку [а, 6], в .расчет не принимаются. При
достаточно больших N среднее арифметическое весьма
близко к математическому ожиданию случайной
величины. Поэтому в качестве приближенного значения для
интеграла A.13) может быть взято среднее
арифметическое
sg
A.14)
Как осуществить это практически? Пусть перед нами
стоит задача вычислить
/= } h(x)dx. A.15)
a
Выберем некоторый закон распределения f(x), для
которого имеется удобный способ получения случайных
чисел и областью определения которого является
интервал (а, Ь). Преобразуем подынтегральную функцию
следующим образом:
U'"S-f(x)dx. A.16)
а
Если теперь обозначить
18

то интеграл примет вид
ь
1= lg*(x)f{x)d(x) A.18)
a
и может быть вычислен при помощи метода
статистических испытаний.
В частном случае, если а и Ь конечны или их можно
считать конечными приближенно, в качестве f(x)
целесообразно выбрать равномерный закон распределения.
Как известно (см. [5]), плотность вероятности
равномерного закона распределения в интервале (а, Ь) равна:
Подставим в интеграл A.16) значение f(x) из
формулы A.19)
ъ
^-dx A.20)
и рассмотрим процедуру вычисления, описанную выше,
для общего случая в связи с равенством A.14).
Из множества равномерно распределенных
случайных чисел выбирается х*. Для каждого значения х\
вычисляется h(xi), затем вычисляется среднее значение
1г{Хг)= — 2 h(Xi) A.21)
функции h(x) в интервале (а, Ь). Таким образом,
величина интеграла A.20) может быть представлена
следующей формулой
I&(b-a)h(Xi). A.22)
На рис. 5 видно, что площадь под кривой
заменяется площадью прямоугольника с основанием (Ь—а) и
высотой, равной среднему значению функции h(Xi) на
отрезке (а, Ь).
2* 19

Рис. 5.
Рассмотренный частный случай находит широкое
применение для вычисления интегралов методом
статистического моделирования в силу того, что границы
области определения легко могут быть приведены к
пределам интегрирования (а, Ь).
3. Кратные интегралы
Рассмотрим приемы вычисления интеграла методом
статистических испытаний с точки зрения
целесообразности их применения для вычисления многократных
интегралов.
Методику подхода к вычислению многократных
интегралов проследим на примере интеграла
-dxn A.23)
по ограниченной замкнутой области Q, заключенной
внутри м-мерного 'параллелепипеда
i; /=1,2, . ..,л.
A.24)
(L25)
Заменой переменных
хг-=сц+(Ьг-аг
приведем интеграл A.23) к виду
где
20

'*= I у- )f[a\+(bi-ai) уьаг+ (Ь2-а2)у2 .
(О
y2-"dyn. A.26)
Область интегрирования целиком заключена внутри
n-мерного единичного куба.
Для вычисления интеграла применим уже
изложенный метод (см. стр. 12). Изменением масштаба по оси
z = f(xbx2, ..., хп) преобразуем вычисляемый интеграл
7 A.27)
где L — минимальное, а М — максимальное значение F
по области со.
A.28)
О^г^ 1,
где
2=/(хь^2,-^п). A.29)
Будем рассматривать вычисляемый интеграл как
некоторый объем в n-мерном пространстве:
/ = § I ••• J • dxidx2—dxndz. A.30)
у
Координаты случайной точки в этом пространстве
имеют вид совокупности п+\ случайных чисел
#=(#1 Rb-#n,U), A.31);
равномерно распределенных в интервале @,1).
Проведем jV испытаний, состоящих в заполнении /г-мерного
единичного куба равномерно распределенными
случайными точками. Будем проверять принадлежность этих
точек объему V посредством неравенства
U <f(/?i,/?2,-,/?n), A-32)
аналогичного неравенству A.2). По окончании
эксперимента получим искомый объем
_ т
ж
~ N
21

Вычисляя интеграл A.23) как среднее значение
функции, поступаем следующим образом. Выбираем
последовательность п чисел
равномерно распределенных в интервале @,1), причем
точки Rij)t не принадлежащие области о, в расчет не
принимаются. Получив N точек /?(j), принадлежащих
области (о, вычисляем среднее значение
/cp-isM/tf», R?---Rl?) A.34)
функции / по области со. Затем определяем значение
интеграла A.23)
/ =5cofCp. A.35)
где S<a — объем области а).
4. Требуемое количество операций
Рассматривая различные способы вычисления
интеграла методом статистических испытаний, мы ввели
N — число испытаний, оговаривая каждый раз, что N
достаточно велико. Какую же величину N можно
считать достаточной для вычисления интеграла с заданной
точностью? Будем говорить, что равенство а^а имеет
точность е с достоверностью а, если вероятность
/>[(|a-oj)<e]=a. A.36)
Проследим это на примере второго способа
вычисления интеграла на ЭВМ (см. стр. 16). В качестве
приближенного значения / интеграла / мы брали величину
!jt. Если рассматривать /= jj- как случайную величину,
можно вычислить математическое ожидание /
AfG)=Af(-Wp С-37)
22

и дисперсию
P(I)= P Р '» A.38)
тогда средняя квадратичная ошибка
./.
Pi\-P)
Л О-39)
Свяжем величины а и 8 с числом испытаний N при
помощи неравенства Чебышева (см. [5]):
^[(|/-/|)<е]^1-4-' A.40)
Заменим левую часть неравенства A.40)
выражением из формулы A.36)
4
о
подставляя вместо а2 его значения, получаем
РA—Р)
и, решая неравенство относительно N, вычислим
A—р)р
^^-7Г-4т-- A-43)
A—а)е2 v _ '
Надо сказать, что формула A.43), полученная на
основании неравенства Чебышева, дает только верхнюю
грань величины N.
Более точно количество испытаний N можно
подсчитать, учитывая, что в рассматриваемом примере
(см. стр. 16) величина Р^-^ асимптотически (при
N—>-оо) подчиняется нормальному закону
распределения. Основываясь на этом факте, перепишем формулу
A.36)
Р\ lEzZL<t lea. (L44)
u a j
Величина ta выбирается из таблиц нормального
распределения для заданного а. Сопоставив равенства
A.40) и A.44), находим
23

A.45)
или
отсюда
NssI^lElt». A.47)
Этой более точной формулой мы будем
пользоваться в дальнейшем для оценки количества испытаний N.
Теперь рассмотрим приемы вычисления многомерных
интегралов с точки зрения сравнения этих приемов с
обычными кубатурными формулами. Выясним
трудоемкость вычисления кратных интегралов методом
статистических испытаний. Полное число операций
k = Nv . A.48)
где v — число операций, затрачиваемое на вычисление
одного значения подынтегральной функции
f(R?,Rl->RZ) A.49)
(предположим, что остальные вычисления
незначительны).
Число операций также зависит от того, какую
точность мы хотим получить, >и. следовательно,
k=k{&). A.50)
Обычно используются кубатурные формулы вида
0 0 0
=Alf(Ql)+A2f(Q2)+-+Arf(Qr), A.51)
где Qi, Q2, ..., Qr — точки n-мерного куба, для
вычисления которых требуются L=rv операций, где L тоже
зависит от заданной точности результата
L = L(e). A.52)
Очевидно, что сравнивать объемы вычислений
методом статистических испытаний и с помощью кубатур-
24

йых формул необходимо соблюдая одинаковую
точность <8.
Выясним сначала зависимость k(e). На основании
формулы A.47)
1р)р 9
-~--t2. A.53)
Если говорить о «максимальной» ошибке, то для
данного случая можно воспользоБаться правилом «трех
сигм». Тогда ta = 3.
Так как максимум значения р(\—р) достигается
при /? = 0,5 и равен 0,25, можно приблизительно оцепить
величину N:
^-> о-54)
следовательно, k имеет порядок
для интегралов любой кратности.
Для кубатурных формул
L(e)~e"v , A.56)
где q зависит от гладкости функции, an — число
измерений. Таким образом, отношение числа операций
&(е), необходимых при вычислении интеграла методом
статистических испытаний, к числу операций Ь(г) при
счете по кубатурным формулам составляет
^~ 2< A.57)
L(e)
Из формулы A.57) следует, что метод
статистических испытаний имеет преимущества уже при
— >2 A.58)
для малых е.
Для наглядности приведем числовые данные о
требуемом количестве испытаний.
25

Пусть достоверность а = 0,95, точность е принимает
значения 0,05; 0,02 и 0,01. Тогда количество испытаний
N при различных значениях вероятности р
представлено в табл. 1. Для каждого е в левой колонке содержится
N, вычисленное по формуле A.43), а в правой — по
формуле A.53), где для а = 0,95 в соответствии с
таблицами нормального распределения /а=1,96.
Табл ица 1
р ^^\
0,1 0,9
0,2 0,8
0,3 0,7
0,4 0,6
0,5
0,05
A 43)
720
1280
1680
1920
200
A.53)
140
250
330
380
390
8,02
A.43)
4 500
8С00
10 Г,00
12 000
12 500
A53)
9С0
1500
2100
>2;ЗС0
'2 400
0,01
A.43)
18С00
32000
46000
48 000
50 000
A.53)
3600
6 200
8 400
9400
9 800
Из рассмотрения таблицы видно, что при переходе от
*р = 0,1 @,9) к р=0,5 /V возрастает примерно в 3 раза, а
при переходе от е = 0,05 к е = 0,01 возрастает примерно
в 25 раз.
В табл. 2 приводятся сравнительные данные о
количестве операций, необходимых для вычисления
кратного интеграла по кубатурным формулам и методом
статистических испытаний, в зависимости от кратности
интеграла п.
Таблица 2
л
2
4
6
8
10
Кубатурные формулы
4-Юз
8-105
1,2- 10s
1,6-lO1©
2.1012
Метод статистических
испытаний
2-105
4-105
6-105
8-105
106
При расчете этой таблицы условно принято, что
интервал интегрирования разбит на 10 частей по каждой
из п осей координат.
Таким образом, метод статистических испытаний
оказывается уже более выгодным при кратности интеграла,
равной 3 или 4.

Глава II
ПОЛУЧЕНИЕ И ПРЕОБРАЗОВАНИЕ
СЛУЧАЙНЫХ ЧИСЕЛ
5. Случайные числа
Как уже говорилось, для реализации метода
статистических испытаний требуется большое
количество случайных чисел. Так как практическое
использование метода статистических испытаний сопряжено с
применением ЭВМ, целесообразно поставить вопрос о
способах получения на ЭВМ достаточно длинных
последовательностей случайных чисел с заданным законом
распределения.
Распространены два основных принципа получения
случайных чисел. При первом из них случайные числа
вырабатываются специальной электронной приставкой
(датчиком случайных чисел), устанавливаемой на ЭВМ.
Реализация этого принципа почти не требует
дополнительных операций машины, кроме операции обращения
к датчику.
Второй принцип — алгоритмический — основан на
формировании случайных чисел в самой машине
посредством специальных программ. Недостатком
алгоритмического метода по сравнению с аппаратурным является
дополнительный расход машинного времени, так как в
этом случае машина сама выполняет операции
небольшой электронной приставки.
Программы выработки случайных чисел с
некоторыми законами распределения могут оказаться
достаточно громоздкими. Поэтому случайные числа с заданным
законом распределения обычно получают не непосред-
27

ственно, а путем преобразования случайных чисел,
имеющих некоторое исходное распределение.
К исходному распределению предъявляются
следующие требования: простота получения чисел с помощью
электронных приставок или непосредственно на ЭВМ,
а также удобство преобразования в распределение с
заданным законом распределения.
На практике считается, что равномерный закон
распределения в достаточной степени удовлетворяет этим
требованиям.
При таком подходе к получению случайных чисел
первую часть работы — выработку случайных чисел с
равномерным законом распределения — выгодно
выполнять аппаратурным методом, так как для этого
достаточно единственной электронной приставки, которая
освобождает ЭВМ от наиболее трудоемкой части
вычислений. В некоторых случаях, когда на данной ЭВМ
задачи решаются методом статистических испытаний лишь
эпизодически и строить электронную приставку
нецелесообразно, случайные числа с равномерным законом
распределения можно вырабатывать и на ЭВМ по
специальным программам.
Напомним свойства равномерного распределения.
Непрерывная случайная величина | имеет
равномерное распределение в интервале (a, b), если ее функция
плотности
1 0 х<а; х^Ь,
а функция распределения
0 х<а
х—а
— а^х^Ь B.2)
1 х>Ь
Это распределение и нужно получить на ЭВМ. Но
оказывается, что (получить точно такое распределение
на ЭВМ невозможно в силу хотя бы того
обстоятельства, что ЭВМ, оперирующая с п-разрядными двоичными
числами, может формировать не более чем 2П
различных чисел, а при равномерном распределении предпола-
28

гается, что бесчисленное множество возможных
значений случайной величины | заполняет непрерывно
интервал (а,Ь). Будем использовать вместо непрерывной
совокупности равномерных случайных чисел интервала
[0,1], дискретную совокупность 2П случайных чисел того
же интервала, имеющих вероятности 2~п. Закон
распределения такой совокупности носит название
квазиравномерного закона распределения. Следует отметить,
что при достаточно большом п различие между
равномерным и квазиравномерным распределением можно
считать практически несущественным. На практике уже
при п, равном 20, различие в случайных числах с этими
распределениями не оказывается заметно на точности
решения задач методом статистических испытаний.
Рассмотрим известные принципы получения
последовательности квазиравномерных случайных чисел:
1) генерирование случайных чисел специальной
электронной приставкой путем моделирования
некоторых случайных процессов;
2) 'получение та-к называемых псевдослучайных
чисел с помощью специального алгоритма.
Идея генерирования случайных чисел, подчиненных
квазиравномерному закону распределения на отрезке
[0, 1], предполагает следующее.
Для получения /г-знач.ного двоичного случайного
числа моделируется последовательность независимых
случайных величин zi} принимающих значения 0 или 1
с равной вероятностью. Полученная последовательность
нулей и единиц представляет собой случайное двоичное
число, квазиравномерно распределенное на отрезке
[0,1].
Аппаратурные методы получения
квазиравномерного распределения различаются только способами
получения последовательности независимых случайных
величин Z{.
Один из возможных способов основан на подсчете
количества радиоактивных частиц за определенный
промежуток времени At. Если число частиц за время
Д^ четное, zi присваивается значение 1, если
нечетное — 0.
Другой способ 'получения случайных величин z{
использует шумовой эффект электронной лампы.
Собственный шум электронной лампы выражается некото-
29

рым выходным напряжением u(t), которое является
случайной функцией. Фиксируя значения этого напряжения
в определенные моменты времени /,-, получаем
последовательность независимых случайных величин u(U).
Пользуясь значениями u(ti), можем получать z\
следующим образом:
.= f 0; u(ti)^a B.3)
"l I 1; u(U)>a
где а выбирается с таким расчетом, чтобы P(zi=\) =
Кроме недостатков аппаратурного принципа
получения случайных чисел, изложенные здесь способы имеют
свои специфические недостатки. Способ первый
предполагает наличие незатухающего или малозатухающего
источника радиоактивных частиц, а второй —
стационарного режима подачи напряжения, позволяющего
сохранять постоянной величину Р(и(и)^а), в противном
случае снижается качество случайных чисел.
Кроме того, аппаратурный принцип не позволяет
использовать для контроля работы программы так
называемый метод двойного просчета, сущность которого
заключается в двойной реализации разбитого на части
алгоритма и последующего сравнения результатов.
Очевидно, что при повторном генерировании не
удается получить те же случайные числа.
Псевдослучайными называются числа,
сформированные на ЭВМ с помощью специальных программ
рекуррентным способом: каждое случайное число получается
из предыдущего с помощью определенных
преобразований. Сформированная последовательность чисел
должна хорошо приближаться к последовательности
случайных чисел с заданным законом распределения.
Рассмотрим простейший (прием получения пссвдосау-
чайных квазиравномерных чисел рекуррентным способом.
Пусть имеется некоторое n-разрядное двоичное число
в интервале [0, 1]. Возведем его в квадрат. Очевидно, что
при этом мы 'получим уже 2/г-разрядное число.
Выделим средние п разрядов этого числа. Полученное
таким образом новое n-разрядное число опять возведем
в квадрат и т. д. Образованная таким образом
последовательность называется псевдослучайной, так как слу-
30

чайность в теоретико-вероятностном смысле здесь Me-
ста не имеет. Однако статистическая проверка
показывает, что полученная этим способом последовательность
чисел близка к последовательности случайных чисел с
равномерным распределением.
Можно привести многочисленные примеры
аналогичных приемов получения псевдослучайных чисел, однако
практическое применение получили лишь те немногие
способы, которые приводят к малым затратам
машинного времени на получение каждого числа.
Следует отметить, что программы для получения
псевдослучайных чисел не всегда работают достаточно
устойчиво и надежно. Во-первых, рекуррентный процесс
может оборваться (выродиться), например,
получившийся при очередном просчете нуль повлечет за собой
нулевую последовательность. Кроме того,
последовательность случайных чисел может оказаться периодической
(ведь число различных чисел не более чем 2V). При
использовании большого количества случайных чисел,
полученных рекуррентным способом, по упомянутым
здесь причинам величины статистических характеристик
их могут существенно исказиться. Для ослабления
влияния возмущающих факторов существует несколько
способов, например, применение периодически
работающих спаренных программ с различными алгоритмами,
где случайное число, полученное с помощью одного из
алгоритмов, является исходным для другого; если же
обращаться к парной программе случайным образом,
то это дает еще больший эффект. Полученные таким
образом случайные числа требуют проверки .на
«случайность», на «близость» к заданному закону
распределения, на отсутствие корреляции и т. д., после чего
они пригодны к использованию.
6. Получение случайных чисел с заданным
законом распределения
Случайные числа (квазиравномерные «и
псевдослучайные с равномерным законом распределения), хотя и
являются равномерными лишь приближенно, могут
быть использованы в качестве исходного материала для
получения любых вероятностных объектов.
31

Такими вероятностными объектами в первую
очередь являются случайные события, наступающие с
заданной вероятностью, случайные величины с заданным
законом распределения и некоторые виды случайных
векторов и процессов.
Посмотрим, как можно моделировать с помощью
случайных равномерно распределенных чисел случайные
события, наступающие с заданной вероятностью. Эту
процедуру называют еще «реализацией жребия». Пусть
событие А наступает с вероятностью р, тогда процедура
моделирования этого события с помощью равномерно
распределенных в интервале @,1) случайных чисел
выглядит следующим образом:
1) выбирается очередное случайное число ?*;
2) проверкой неравенства
Ь^Р B.4)
устанавливается принадлежность этого числа отрезку
[0, р].
Если число | удовлетворяет неравенству B.4),
говорят, что событие А наступило, в противном случае —
не наступило.
Аналогично выглядит процедура моделирования на
ЭВМ дискретной случайной величины с заданным
законом распределения.
Пусть случайная величина ? принимает возможные
значения Z\, z% ..., zn с вероятностями р\, р2> -•, Рп.
Очевидно, что значение Z\ будет принято случайной
величиной ? в том случае, когда выполняется
неравенство, аналогичное неравенству B.4)
1<<р; B.5)
(наступает событие, состоящее в том, что ?=--:?i);
значение 22 — когда
Pi<li<P\+P2 B.6)
(наступает событие, состоящее в том, что ? = 22);
значение ^з — когда
Р\+ р2^Ь<Р\ + Р2 + Рз B.7)
(наступает событие, состоящее в том, что 'Q^2*) и т. д.
Другими словами, пусть

tr= 2 pi • B.8)
Тогда, если
lmr-1 <KU B.9)
наступает событие, состоящее в том, что Z> = zm.
Процедура реализации этого способа моделирования
дискретной случайной величины на ЭВМ сводится к
следующему. Вырабатываем случайные числа |* с
равномерным распределением в интервале @,1).
Очередное |г сравниваем с 1\\ если неравенство B.5) выполнено,
считаем, что ?=Z\; в противном случае переходим к ^.
Сравниваем & с W, если неравенство B.6) выполнено,
считаем, что ? = 22; в противном случае переходим к /3
и т. д. до тех пор, пока одно из неравенств вида B.9)
окажется выполненным. Эта процедура всегда рано или
поздно приводит к цели, так как событие, состоящее в
в том, что случайная величина | принимает какое-аш-
будь из своих значений zit является достоверным.
Перейдем к рассмотрению метода моделирования
непрерывных случайных величин. Пусть по-прежнему в
нашем распоряжении имеются случайные числа & с
равномерным распределением в интервале @,1).
Требуется получить случайные числа у и являющиеся
возможными значениями случайной величины ц с законом
распределения, заданным функцией плотности / (у).
Можно доказать, (см., например [9]), что случайная
величина г), являющаяся решением уравнения
B.10)
— 00
имеет распределение /л (у), если случайная величина I
распределена равномерно в интервале @,1).
Соотношением B.10) можно воспользоваться для
получения случайных чисел с заданным законом
распределения.
Методику преобразования случайных чисел
поясним па примерах.
Пример 1. Предположим, что нам необходимо
получить случайные числа с показательным
распределением
B.11)
3 Заказ 3012. 33

Используем соотношение B.10)
e-?ydy = li • B.12)
о
Интеграл B.12) берется в конечном виде, поэтому
1-е = ?,. B.13)
Решим уравнение B.13) относительно //*:
Уь=-\ In (I—gi). B.14)
Соотношение B.14) полностью решает поставленную
задачу. Заметим, что случайное число |'=1—? имеет
также равномерное распределение в интервале @,1).
Поэтому вместо B.14) обычно пользуются
соотношением
tji^lnh • B.15)
Подставляя в правую часть соотношения B.15)
последовательно случайные числа ii, мы получим
последовательность чисел t/i с показательным законом
распределения.
Пример 2. Случайная величина г\ имеет функцию
плотности
M*/)=*(i--^)' (°<у<т) Bл6)
Располагая случайными числами |г-, имеющими
равномерное распределение в интервале (ОД), требуется
получить случайные числа с законом распределения,
который выражается функцией плотности B.16).
Соотношение B.10) имеет вид:
<2Л7>
о
После вычисления интеграла
Я [Ui~~y2i ) =!;• B-18)
34

Решив уравнение B.18) относительно у и получим:
у.= L (i-yT^) B.19)
или, с учетом замечания, относящегося к соотношению
B.14):
0i=4-O-V&)" B.20)
Можно привести и другие примеры использования
соотношения B.10). Однако приходится признать, что
практически изложенная методика имеет весьма
ограниченную сферу применения. Это объясняется
следующими двумя обстоятельствами.
Во-первых, для многих законов распределения,
встречающихся в практических задачах, интеграл B.10)
в конечном виде не берется. Например, для
нормального распределения (закона Гаусса)
1 <*-*?.
а* B.21)
1'2л а
соотношение B.10) приводит к интегралу
—rr- А е
\'2ла
который можно вычислить только численными
методами. Это приводит к недопустимо большим затратам
машинных операций на преобразование случайных чисел.
Во-вторых, даже для тех случаев, когда
соответствующие интегралы берутся в конечном виде,
получаются формулы, например B.15), B.20) и др., весьма
неудобные для расчета на ЭВМ. Причина в том, что
вычисление логарифмов, корней и других элементарных
функций на ЭВМ выполняется при помощи
стандартных программ, состоящих из многих исходных
операций машины (сложение, умножение и т. д.). Учитывая,
что применение метода статистических испытаний
требует большого количества случайных чисел,
естественно, возникает опасность существенных затрат машин-
3* 35

ного времени. Поэтому на практике обычно пользуются
приближенными методами преобразования случайных
чисел, обеспечивающими достаточную экономию
операций ЭВМ.
7. Приближенные способы преобразования
случайных чисел
Рассмотренная процедура моделирования случайной
величины с заданным законом распределения с
помощью равномерно распределенной случайной величины,
имеет существенные недостатки, препятствующие
широкому применению этого метода на практике. Это
обстоятельство послужило толчком для создания
различных способов приближенного моделирования.
Существуют способы, пригодные только для
моделирования случайных величин с конкретными законами
распределения, существуют также универсальные способы,
с помощью которых возможно моделировать законы
распределения любого вида.
Приведем один из таких универсальных способов.
Пусть закон распределения случайной величины т|,
предложенной нам для моделирования, задай функцией
плотности /г, (у), возможные значения которой лежат в
интервале (а, Ь). Если это интервал с бесконечными
границами, целесообразно перейти к усеченному
распределению. Представим / (у) на участке (а, Ь) в виде
кусочно-постоянной функции (см. рис. 6), т. е. разобьем
(а, Ь) на п интервалов и будем считать fn {у) на каж-
О о
Рис. 6.
36

дом интервале постоянной; тогда случайную величину
г] можно представить в виде
ri=ah + lb B.23)
т. е. на каждом участке (аи, ал+i) величина |/г считается
распределенной равномерно. Чтобы аппроксимировать
/л (у) наиболее удобным способом, целесообразно
разбить (af b) на интервалы таким образом, чтобы
вероятность попадания случайной величины г) в любой
интервал (ад, a/i+i) была постоянной, т. е. не зависела от
номера интервала k.
Для вычисления аи пользуются следующим
соотношением:
п
B.24)
где п — количество интервалов (обычно принимается
равным 2т, где т — целое положительное число).
Процедура моделирования предполагает следующее:
1) выбирается случайное равномерно
распределенное число |г-;
2) с помощью |г случайным образом выбирается
интервал (aft, fl/i+0;
3) берется следующее равномерно распределенное
число |i+i и масштабируется с целью приведения его к
интервалу (ak, flfc+i), т. е. |,-+i становится случайной
величиной, равномерно распределенной в интервале {аи,
)
+)
Случайное число у% с требуемым законом
распределения вычисляем по формуле
yi = ah+ l2+i {ak+i - ah). B.25)
Интересно несколько подробнее рассмотреть процесс
выборки интервала (аи, fl/i+i) с помощью |^. Как
известно, случайное число & получается на машине в виде
последовательности нулей и единиц. Если число
интервалов п равно 2т, то количество всевозможных
комбинаций нулей и единиц m-разрядного двоичного числа
даст нам количество интервалов п, т. е. каждому
интервалу (аи, CLk+\) можно поставить в соответствие одну и
только одну комбинацию нулей и единиц т-разрядного
37

двоичного числа. Это обстоятельство сильно упрощает
методику выбора интервала (аи, Яй+О-
Каждый интервал кодируется m-разрядным
двоичным числом, и результаты заносятся в специальную
таблицу.
В эту же таблицу заносятся значения пи—ah~\ —
коэффициента масштабирования.
Таким образом, выделив т разрядов ?*, мы сразу
определяем номер интервала (а&, au+i).
Для реализации изложенного метода приближенного
моделирования случайных величин на ЭВМ требуется
небольшое количество операций; кроме того, количество
операций не зависит от точности аппроксимации (т. е.
от количества интервалов п). Точность аппроксимации
влияет только на размеры участка памяти, куда
помещается таблица закодированных значений а/?. Этим
способом преобразования случайных чисел широко
пользуются в практике статистического моделирования.
Небезынтересны также способы получения
случайных чисел с заданным законом распределения,
основанные на использовании предельных теорем теории
вероятностей. В качестве примера рассмотрим способ
получения случайных чисел с нормальным законом
распределения, весьма часто встречающийся при решении
практических задач.
В силу центральной предельной теоремы теории
вероятностей сумма большого числа одинаково
распределенных независимых случайных величин при весьма
общих условиях имеет приближенно нормальное
распределение.
Пусть f ь &». • • * in случайные величины, имеющие
равномерное распределение в интервале @,1). Тогда
- -+1п B.26)
оказывается случайной величиной с распределением,
близким к нормальному, при больших п.
Как известно, математическое ожидание (среднее
значение) случайной величины |<
B.27)
а среднее квадратическое отклонение
38

——• B-28)
2УЗ
Поэтому математическое ожидание а суммы г)
B.29)
П
а —-,
а среднее квадратическое отклонение
у7
Процедура формирования случайных чисел с
нормальным распределением, имеющим а=0 и сг=1,
сводится к следующему:
1) выбираются п последовательных случайных чисел
? ? 1
+ +
2) ВЫЧИСЛЯеТСЯ Сумма Г) = ?г + ?г+1+ ••• +^г+л',
3) определяется случайное число
п
Ч —
2 У3Bт|-/г)
т|*= —= _—> B.31)
2](Г
имеющее приближенно нормальное распределение со
средним значением, равным 0, и дисперсией —1. Как
показывает опыт, для решения практических задач
можно пользоваться я = 5-т-6.
8. Моделирование случайных векторов
При решении задач методом статистических
испытаний нередко возникает необходимость в формировании
возможных значений или, как иногда говорят,
реализаций случайных векторов.
Случайный вектор можно задать проекциями на оси
координат, причем эти проекции являются случайными
39

величинами, описываемыми совместным законом
распределения. В простейшем случае, когда
рассматривается случайный вектор на плоскости XOY, необходимо
задать совместный закон распределения его проекций
? и ц на оси X и У соответственно.
Предположим сначала, что двумерная случайная
величина (|, г\) является дискретной и ее составляющая
| принимает возможные значения Х\, Х2,..., хп> а
составляющая ц — значения у и \)ч, • . •, \}п> причем каждой
паре (xi, yj) соответствует вероятность рц. При этих
предположениях молено найти частное распределение
случайной величины |, а именно: каждому возможному
з-начению х-ь случайной величины § будет
соответствовать вероятность
/?г= 2 Piy B.32)
Теперь по правилам, рассмотренным в § 6 (см. B.8),
B.9) и др.), можно определить конкретное значение х\
случайной величины | в соответствии с распределением
вероятностей B.32). Пусть это будет хц. Тогда из всех
значений ри выберем совокупность
Pi^ Pip • ¦ • . Pivn> B.33)
которая описывает условное распределение случайной
величины г] при условии, что |=хг\ Затем по тем же
правилам определим конкретное значение у\ (пусть оно
равно у и) случайной величины г| в соответствии с
распределением вероятностей B.33). Полученная пара
(хц\ У и) и будет первой реализацией моделируемого
случайного вектора. Далее аналогичным способом
определяем возможное значение Х\2 в соответствии с
распределением B.32), выбираем совокупность
и находим у 12 в соответствии с распределением
вероятностей B.34) и т. д.
Процедура моделирования не претерпевает
принципиальных изменений и.в том случае, когда речь идет о
моделировании непрерывного случайного вектора. В
этом случае двумерная случайная величина (|, ц)
описывается совместной функцией плотности f (х, у). Част-
40

пая функция плотности случайной величины ? может
быть определена в соответствии с соотношением,
аналогичным B.32), а именно:
B.35)
—оо
Имея функцию плотности f% (x), можно «найти
случайное число %г (пусть это будет х\) по правилам,
рассмотренным в § 7. Затем определяется условное
распределение случайной величины г) при условии, что ?=Х\.
h №) B.36)
В соответствии с функцией плотности B.36) можно
определить случайное число yi (пусть это будет у{). Тогда
пара (х, у) и является искомой реализацией вектора
(I Л).
Рассматриваемый метод формирования реализаций
случайных векторов в принципе обобщается на случай
пространства любого числа измерений п (>на многомерные
случайные величины). Однако при больших п объем
вычислений существенно возрастает и может служить
серьезным препятствием для практической работы.
Процедура формирования реализаций случайного
вектора значительно упрощается и становится мене&
громоздкой в том случае, когда многомерная случайная
величина задается в рамках корреляционной теории
(при помощи корреляционной матрицы).
Остановимся кратко на трехмерном случае. Пусть
требуется сформировать реализации трехмерного
случайного вектора (?, rj, ?), имеющего нормальное
распределение с математическими ожиданиями
Af(g)=ai; Af(n)=fl2;
а3 B.37)
и корреляционной матрицей
Kii К" 12 KiZ
Кi2 Кгг Кгз
С2 3
А13 Л 32 А 33
41

Здесь /Си, /С22, /Сзз — дисперсии случайных величин
?, Л и ? соответственно, Ki2 = K2i, Ki3= K31 и К2з=Кз2—
корреляционные моменты | и л» Е и ?, р?
соответственно.
Будем предполагать, что в нашем распоряжении
имеются случайные числа v*, имеющие одномерное
нормальное распределение с математическим ожиданием т
и дисперсией а2. Способы формирования случайных
чисел, обладающих такими свойствами, мы уже
рассматривали.
Выберем три числа, скажем vi, V2, v3, и преобразуем
их так, чтобы они имели характеристики B.37) и B.38).
Искомые составляющие случайного вектора (|, л>
?) обозначим х, у, z к представим в виде:
- m)+c22{v2-rn)+a2\
—m) + с23(\«-т) +c33(v3-m) +a3, B.39)
где Cij — пока неизвестные нам коэффициенты. Для
вычисления коэффициентов сц воспользуемся
элементами корреляционной матрицы. По определению
но МA)=аи
M[(vi-m)]2 = a2;
поэтому
Ки = с*цо\ B.40)
Аналогично, поскольку М [(v» — т) (vj — т)] = 0 при
[ф] (случайные величины v\, L2, V3 независимы между
собой):
/C 2
Кзз = Ci3 a2+с22з о2 + cl a2. B.41)
42

Соотношения B.40) и 2.41) 'представляют собой
систему уравнений относительно коэффициентов сц. Ре
шив эту систему уравнений, найдем:
а
Kit
V /С11/С22—/
-К.
?'22 =
с13=—-г=; B.42)
а у/Си
У /Cii'/Сгз—/Ci2*Ai3
v
ау/Сп
Располагая коэффициентами Cij B.42), легко три
последовательных независимых случайных числа v* пре-
преобразовать в составляющие случайного вектора (Хг\
Уг\ Zi) вида B.39).
На этом мы закончим рассмотрение приемов
моделирования элементарных вероятностных схем при
помощи случайных чисел. Аналогичные приемы могут быть
разработаны и для других случаев, встречающихся при
решении практических задач.

Глава III
ПРИМЕНЕНИЕ МЕТОДА СТАТИСТИЧЕСКИХ
ИСПЫТАНИЙ ДЛЯ МОДЕЛИРОВАНИЯ
СЛОЖНЫХ СИСТЕМ
9. Математическая модель
Рассмотренный выше метод статистических
испытаний может быть использован для моделирования
не только простейших вероятностных схем, ной процессов
функционирования различных реальных систем,
встречающихся в технологии, организации производства,
автоматическом управлении, экономике, планировании,
учете и обработке информации.
Моделирование процесса возможно лишь в том
случае, если для него построено четкое формальное
описание, учитывающее основные закономерности процесса и
действующие факторы.
Процесс представляет собой последовательную
смену состояний системы во времени. Каждое мгновенное
состояние системы описывается набором чисел zu
22,..., 2n-i, zn, которые достаточно полно отображают
основные свойства системы в данный момент. Очевидно,
что в общем случае значения г{ зависят от времени и,
следовательно, могут быть выражены как функции
времени Z\ (/), z2 (/),..., zn (t). Эти величины будем
называть в дальнейшем характеристиками состояний
процесса. Значения zu 22,..., zn можно интерпретировать
как координаты точки в /г-мерном фазовом
пространстве. Если рассматривать фазовую координату z{ (t) как
функцию времени /, нетрудно представить фазовую
44

траекторию .как вектор-функцию z(t) с составляющими
по осям координат г{ (/), z2 (/),... ,гп (/).
Математическая модель реального процесса есть
некоторый математический объект, поставленный в
соответствие данному физическому процессу. В дальнейшем
под математической моделью будем понимать
совокупность соотношений, связывающих характеристики
состоянии процесса с параметрами системы, исходной
информацией и начальными условиями. Это, однако, не
означает, что математическая модель состоит только из
соотношений, выражающих характеристики состояний
процесса как явные функции параметров, начальных
условий и исходной информации. В общем случае этого
может и не быть. Сущность математической модели
заключается в том, что при совместном рассмотрении
составляющих ее соотношений характеристики состояний
процесса однозначно определяются как функции
упомянутых выше аргументов. Однако однозначность (в
полном смысле слова) определения характеристик
состояний имеет место только для вполне детерминированных
процессов, при исследовании которых не учитываются
случайные факторы. На практике же чаще всего
приходится изучать случайные процессы, характеристики
состояний которых описываются случайными функциями
времени.
Характеристики состояний процесса могут быть
случайными функциями в силу различных обстоятельств:
вполне детерминированные системы могут иметь
случайные начальные условия или случайные параметры;
при вполне детерминированных параметрах и
начальных условиях сами процессы также могут быть
случайными. Последнее обстоятельство чаще всего
объясняется воздействием на элементы системы случайных
возмущений, возникающих как внутри, так и вне данной
системы. Наиболее распространенной на практике
является ситуация, когда сам процесс, начальные условия
и параметры описываются случайными функциями.
Рассматривая этот наиболее общий случай, будем
говорить, что с помощью математической модели
однозначно определяется распределение вероятностей для
характеристик состояний процесса, если задано
распределение вероятностей для начальных условий, параметров
системы и возмущений, действующих на ее элементы.
45

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

значительно шире, чем класс уравнений, доступный
аналитическому исследованию. Однако исследование
процесса с помощью численных методов оказывается, как
правило, менее полным по сравнению с аналитическим.
В качестве результатов при исследовании процессов
численными методами обычно получают таблицы з-ча-
чоний искомых величин для конечного набора значений
параметров системы. Применение средств
вычислительной техники (в том числе и быстродействующих ЭВМ)
при данном способе использования математической
модели ограничивается лишь автоматизацией
вычислений — автоматическим воспроизведением выбранного
численного метода.
В этом отношении моделирование процесса на ЭВМ
принципиально отличается от исследования его
аналитическими или численными методами.
При моделировании процесса па ЭВМ имеет место
воспроизведение происходящих явлений с сохранением
их логической структуры и расположения во времени.
Для моделирования процесса на ЭВМ необходимо,
преобразовать его математическую модель в так
(Называемый моделирующий алгоритм. Реализация
моделирующего алгоритма на ЭВМ является как бы имитацией
явлений исследуемого процесса с учетом их
взаимодействия. Легко заметить некоторую аналогию между
моделированием процесса на ЭВМ и экспериментальным
исследованием процесса (в натуре, на макете,
лабораторной установке и т. д.). В том и в другом случае
последовательно воспроизводятся состояния процесса (при
экспериментальном исследовании — физически, при
моделировании на ЭВМ — путем вычисления координат
zi, Z2, ...,zn), наблюдение и фиксация которых в
нужные моменты времени позволяют получить сведения,
необходимые для исследования процесса.
10. Этапы формализации реальных сложных систем
Рассмотрим в общих чертах методику формализации
сложных систем.
Очевидно, что для сложных систем составление
математической модели непосредственно по результатам
47

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

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

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

Остановимся подробнее на методике изображения
математической модели в операторной форме. Как
следует из названия, операторная схема состоит из
операторов. Существенны два класса операторов:
арифметические операторы и логические операторы.
Арифметический оператор представляет собой в
определенном смысле замкнутую группу элементарных
вычислительных актов и описывает обычно
совокупность каких-нибудь соотношений (например,
вычисление по формуле). Будем обозначать арифметические
операторы заглавными буквами с индексами,
указывающими номер оператора (например, А25). В блок-
схеме моделирующего алгоритма (графической
форме представления операторной схемы) арифметический
оператор обозначается прямоугольником с записанным
внутри него реализуемым соотношением. Передача
управления от арифметического оператора обозначается
индексом справа вверху (например, А25 означает, что
от оператора № 25 управление передается оператору
№ 31).
Логический оператор описывает проверку условий
вида а^.Ь и в зависимости от исхода этой проверки
формулирует признак (со=1, если условие выполнено, и
со = 0, если оно не выполнено). Основное отличие
логического оператора от арифметического заключается в
том, что логический оператор предопределяет
разветвление вычислительного процесса: управление
от'логического оператора может быть передано одному из двух
различных операторов в зависимости от значения
признака со. Будем обозначать логический оператор буквой
Р с соответствующим номером, например Pis. Передача
управления от логического оператора обозначается
следующим образом. Справа от символа логического
оператора ставится стрелка с номером того оператора,
которому передается управление. По стрелке,
направленной вверх, управление передается в том случае,
когда условие, проверяемое данным логическим
оператором, выполнено. Например, Pit6 означает, что если
условие, проверяемое оператором Pi2> выполнено, то
управление передается оператору № 6. Аналогично по
стрелке, направленной вниз, управление передается в
том случае, когда условие, проверяемое данным
логическим оператором, не выполнено. Например, P^ls

означает, что если условие, проверяемое оператором
Pi2> не выполнено, то управление передается оператору
№ 8. Таким образом, в операторной схеме логический
оператор Р}2 выглядит как Рщв ,
В блок-схемах логические операторы изображаются
в виде кружков, внутри которых записываются
проверяемые условия, а выходящие стрелки имеют индексы
О и 1. По стрелке с индексом 1 управление передается
в том случае, когда проверяемое условие выполнено, а
по стрелке с индексом 0, если оно не выполнено.
Как для арифметических, так и для логических
операторов будем придерживаться правила: изображение
v передачи управления от данного оператора оператору,
непосредственно следующему за ним, опускается.
Для обозначения передачи управления данному
оператору от других операторов номера последних
записываются в виде индексов слева вверху от символа
данного оператора, например 15-21А22 означает, что
оператору А22 управление передается от операторов
№ 15 и 21. Передача управления данному оператору от
предыдущего изображается лишь в том случае, когда
данному оператору передается управление от
нескольких операторов.
Арифметические операторы составляют весьма
обширный класс операторов, который подразделяется на
несколько подклассов, соответствующих содержанию
операций, описываемых оператором. К наиболее часто
употребляемым можно отнести операторы следующих
подклассов: вычислительные (обозначаются А),
операторы формирования реализаций случайных процессов
(Ф), формирования неслучайных величин (F), счетчики
(К), а также оператор (Я), обозначающий конец
вычислений и выдачу результатов.
Для иллюстрации обозначений рассмотрим пример
операторной схемы алгоритма и его блок-схемы. Пусть
речь идет о решении квадратного уравнения вида
Как известно, корни этого уравнения выражаются
формулой
52

Для представления в операторной форме алгоритма
й
определения корней
Р
Ai — ;вычисление —к
введем следующие операторы:
Р2
— вычисление D= ~т-\ —
Рз— проверка условия ^;
А4 — определение действительных корней
А5 — определение комплексных корней
Яб — выдача результатов.
Операторную схему рассматриваемого алгоритма
можно записать так:
На рис. 7 представлена
достаточно наглядная блок-
схема этого алгоритма.
Работа алгоритма
сводится к следующему. Опера-
р
тор Ai вычисляет к- и
передает управление
оператору Аг.
Оператор А2 вычисляет
Р2
D= -? q и передает
управление оператору Рз.
Оператор Р3 проверяет условие
D^O. Если это условие
выполнено, оператор Рз
вырабатывает со = 1, в противном
случае о = 0.
Предположим, что D^O и (о= 1.
Тогда от оператора Р3 по
стрелке с индексом 1
управление передается оператору
А4. Оператор А4 вычисляет
действительные корни урав-
р
Вычисление - -у
Вычисление D- а
/у
>,о\
в
Выдача хп
Рис. 7.
53

р —
нения #i,2 = 2~±У|?>| ипередает управление
оператору Яб, который выдает готовый результат.
Возвратимся теперь к оператору Р3 и предположим, что Z)<0.
Тогда от оператора Р3 по стрелке с индексом 0
управление передается оператору А5. Оператор А5
определяет комплексные корни уравнения х\г2 = —2~±/y|D| и
передает управление оператору Яб для выдачи
готового результата.
В дальнейшем мы познакомимся с более сложными
алгоритмами, представляемыми при помощи
операторных схем и соответствующих им блок-схем.
12. Метод статистического моделирования как
аппарат анализа сложных систем
Моделирование процесса функционирования
сложной системы на ЭВМ позволяет провести ряд
количественных и качественных исследований системы,
представляющих значительный интерес 'как с теоретической,
так и с практической точек зрения. По ценности
результаты моделирования в ряде случаев не уступают
результатам эксперимента, вместе с тем, они могут быть
получены при заведомо меньишх затрат времени и средств.
Иногда метод моделирования оказывается более
эффективным инструментом анализа систем, чем
эксперимент. Для примера достаточно привести задачи,
связанные с изучением поведения системы при условиях
или значениях параметров, которые еще не могут быть
созданы в натуре.
Как показывает опыт применения метода
статистического моделирования, полезные результаты,
относящиеся к исследуемой системе, иногда удается получить
уже на этапе построения модели. Необходимость
четкого формального описания процесса
функционирования системы и строгой постановки задач исследования
зачастую являются причиной более глубокого
осмысливания закономерностей, свойственных, системе, и
вскрытия таких подходов к возникающим задачам, которые
были недоступными ранее.

Однако наиболее существенные результаты можно
получить путем реализации модели на ЭВМ. Даже
общий обзор результатов моделирования позволяет
оценить правильность функционирования системы в
различных условиях, вскрыть непредвиденные особенно-
сти в ее поведении.
Богатство сведений, получаемых методом
моделирования об исследуемой системе, можно подтвердить
примером. Если речь идет о моделировании
производственного процесса предприятия, выпускающего штучные
изделия (трубы, автомобили, часы и т. д.), то
величинами, которые фигурируют в качестве результатов
моделирования, обычно считают:
среднее количество готовых изделий, выпускаемых
за смену;
среднюю долю бракованных изделий;
относительное время производительной работы и
простоев каждого станка, агрегата или
технологической линии;
среднее количество нарушений нормального режима
процесса, определяемое для каждой из причин,
вызывающих нарушения;
среднее количество бракованных изделий,
определяемое для каждой из причин, порождающих брак;
среднее время отклонения от графика и нарушения
синхронности процесса для каждой единицы
производственного оборудования;
среднюю продолжительность работы каждого из
видов инструмента;
среднее снижение производительности,
определяемое для каждой из причин, способной снизить
производительность, и т. д.
Аналогичные характеристики можно получить и в
результате моделирования других сложных систем.
Для того чтобы сделать анализ системы методом
моделирования более целенаправленным и четко
сформулировать задачи анализа, обычно выбирают
специальные показатели, характеризующие качество
функционирования системы. В первую очередь речь идет,
как правило, о выборе показателя эффективности
системы— такой числовой характеристики, которая
оценивает степень приспособленности системы к
выполнению стоящих перед ней задач.
55

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

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

Глава IV
СТАТИСТИЧЕСКИЕ МОДЕЛИ СИСТЕМ
МАССОВОГО ОБСЛУЖИВАНИЯ
13. Понятие системы массового обслуживания
Перейдем теперь к анализу методики
моделирования систем некоторого класса, имеющих широкое
практическое применение.
Для каждого типа реальных систем выбирается своя
формализованная схема, способная представить
процесс функционирования системы с достаточной для
практики точностью.
Рассмотрим системы, характеризуемые тем
свойством, что они выполняют некоторые операции над
объектами, поступающими из внешней среды.
В качестве примера рассмотрим работу
бензозаправочной станции. На станцию прибыла машина. Если к
моменту ее появления имеются свободные
бензозаправочные колонки, машина становится на заправку, в
противном случае машина ожидает, пока не
освободится одна из колонок станции.
Процесс работы станции описывается случайными
величинами; поток машин, поступающих на станцию,
является случайным (машины прибывают в случайные
моменты времени), время обслуживания каждой
машины тоже величина случайная (в зависимости от
количества заправляемого бензина) и, наконец, время
ожидания в очереди тоже случайная величина (зависит от
количества машин в очереди). Качество обслуживания
машин на бензозаправочной станции можно
характеризовав средней длиной очереди, средним временем ожи-
58

дания в очереди или средним временем пребывания
машины на станции (с момента прибытия на станцию до
момента окончания заправки).
Аналогично можно описать работу других реальных
систем, например разгрузку пароходов в порту,
обслуживание самолетав взлетно-посадочными полосами,
обслуживание клиентов в парикмахерских и абонентов
на АТС и т. д. Эти системы, различные по своей
физической природе, имеют тем не менее сходную структуру
процессов функционирования. Все они характеризуются
потоком заявок (машин на бензозаправочной станции,
самолетов в аэропорту, клиентов в парикмахерских),
наличием обслуживающих каналов или линий
(бензозаправочных колонок, взлетно-посадочных полос, мастеров
в парикмахерской). Для формального описания систем
такого типа существуют специально разработанные
математические схемы, которые получили название систем
массового обслуживания.
Основными процессами в системах массового
обслуживания являются процесс поступления заявок и
процесс собственно обслуживания заявок каналами
(линиями) системы. Остановимся кратко на
характеристике этих процессов.
Заявки, 'поступающие в систему, образуют поток,
т. е. последовательность событий, специальным образом
расположенных во времени. Если с точки зрения
обслуживания все заявки данного потока равноправны, то
такой поток называют потоком однородных событий.
В этом случае каждое событие характеризуется только
моментом времени t, в который оно 'поступает.
Чтобы описать детерминированный поток,
достаточно задать набор конкретных значений tj.
Для описания случайных потоков однородных
событий задается совместный закон распределения
случайных величин tj; /=1,2, ... k.
Как правило, tj целесообразно заменять величинами
gj, являющимися длинами интервалов времени между
моментами поступления заявок:
D.1)
59

Совокупность случайных величин ?7- задается
совместным законом распределения вероятностей.
Важным классом потоков однородных событий,
имеющих многочисленные .применения, являются
стационарные ординарные 'потоки с ограниченным последствием
(так называемые потоки Пальма). Случайный пото1К
называется ординарным, если вероятность Я(/о, t)
появления двух и более событий за промежуток .времени (t0,
U + t) для любого t является величиной малой по
сравнению с /. Случайный поток называется потоком с
ограниченным последействием, если интервалы ?j между
моментами поступления заявок есть независимые
случайные величины. В этом случае каждый интервал ^ можсг
быть задан своей функцией плотности /j(z). Поток
однородных событий называется стационарным, если его
вероятностный режим не зависит от времени, т. е. все
интервалы ?,-; /= 1, 2, ... между заявками имеют одинаковое
распределение за исключением, быть может, первого
интервала ?ь Для />1 функция плотности ?j, равная
fj(Zj)9 является условной функцией плотности при
условии, что в 'начальный момент интервала ?j поступила
предыдущая заявка. Первый 'интервал ?i отличается от
ц.7 при />1 тем, что относительно поступления или
непоступления заявки в момент f = 0, никаких предположений
не делается. Поэтому f\{z\) представляет собой
безусловную функцию плотности.
Обычно «потоки Пальма задаются функцией плотности
fi(zi) =f(z) при />1. Функция плотности f\(z\) может
быть найдена через f(z).
Прежде чем перейти к рассмотрению
соответствующей формулы, введем понятие интенсивности потока.
Математическое ожидание ?j при />1 равно:
ц= S zf(z)dz. D.2)
0
Здесь ji — средняя длина интервала между
последовательными заявками. Величина
Я-1 D.3)
носит название интенсивности потока и определяет
среднее количество заявок, поступающих за единицу времени.
GO

Для стационарных потоков с ограниченным
последействием имеет место соотношение (формула Пальма),
связывающее функции плотности f\{z\) и f(z):
f(u)du]. D.4)
о
Это соотношение и позволяет получить функцию
плотности /l(Zi) п0 f(Z)-
-Случайный поток однородных событий с
ограниченным последействием может быть потоком без
последействия, если закон распределения оставшейся части
интервала времени между заявками не зависит от того, сколько
этот интервал уже длится. В практике применения
теории массового обслуживания важную роль играет так
называемый простейший поток однородных событий.
Поток называется простейшим, если он является
стационарным, ординарным и потоком без последействия.
Для простейшего потока вероятность Ph(t) наступления
k событий за интервал времени длины / выражается
законом распределения Пауссона.
e-u, D.5)
поэтому простейший поток часто называют пауссонов-
ским.
Функция плотности f(z) случайной величины ?,- при
/>0 для простейшего потока имеет вид показательного
распределения с .параметром Я:
f(z)=%e~te , ,D.6)
где X—интенсивность потока. Другие примеры потоков
Пальма мы рассмотрим в следующем параграфе.
До сих пор мы рассматривали только ординарные
потоки. На практике встречаются случаи обслуживания
групповых заявок, образующих «сгустки» событий. Для
того чтобы описать неординарный поток, необходимо,
кроме момента поступления заявки tj, задать
распределение количества заявок для каждого момента tj.
Рассмотрев методику математического описания
потоков однородных событий, перейдем к изучению
процессов обслуживания заявок.
61

В общем случае система массового обслуживания
состоит из лений (каналов), способных параллельно и
независимо друг от друга обслуживать поступающие в
систему заявки.
Каждая линия может находиться в двух состояниях:
линия свободна или занята. Поступившая в систему за-
-Я'В-ка три наличии свободных линий принимается к
обслуживанию, в противном случае заявка ожидает некоторое
время т, после чего получает отказ и уходит из системы.
В зависимости от величины т системы можно
разделить на тр1и группы:
1. Система с ожиданием (т = оо). Поступившая
заявка терпеливо ожидает возможности быть обслуженной
до тех пор, пока эта возможность ей не представляется
(например, разгрузка пароходов, прибывших в порт
назначения, и др.).
2. Система с отказами (т = 0). При невозможности
быть принятой к обслуживанию немедленно (все линии
заняты) заявка получает отказ и уходит из системы
(например, автоматическая телефонная станция: если
требуемый номер занят—мы слышим короткие гудки).
3. Система смешанного типа (о^т^оо). Это
наиболее распространенный тип систем массового
обслуживания. Поступившая заявка ожидает обслуживания
некоторое время, а затем получает отказ и уходит из системы.
Кроме времени ожидания т, для характеристики
свойств обслуживающей системы обычно вводится
величина т* — длительность обслуживания заяв-ки или время
занятости линии.
Обычно т и т* оказываются случайными и задаются
соответствующими законами распределения, однако
возможны случаи, когда они считаются фиксированными.
Остановимся теперь на возможных видах дисциплины
очереди и порядка обслуживания заявок. Заявки,
ожидающие в очереди, занимают освободившуюся линию в
соответствии со следующими возможными правилами:
1) линию занимает заявка, которая раньше других
поступила в систему; 2) линию занимает заявка, для
которой оставшееся время пребывания в системе наименьшее;
3) заявки принимаются к обслуживанию ,в случайном
порядке; в простейшем случае, если имеется одна
свободная линия и п ожидающих заявок, то вероятность быть
принятой к обслуживанию для -каждой заявки Pt*=—;
62

в более общем случае Pi* вычисляется с учетам
оставшегося времени ожидашия и т. д.
В системах массового обслуживания встречается
обслуживание с преимуществом. Каждой заявке,
поступившей в систему, присваивается некоторый коэффициент
преимущества. С учетом этого обстоятельства к
обслуживанию принимается та заявка, коэффициент
преимущества которой наибольший. Возможны случаи, .когда
запятая линия прекращает обслуживание заявки, если в
систему поступает заявка с большим коэффициентом
преимущества, причем в некоторых случаях недообслу-
женмая заязка получает отказ и уходит из системы, а
для других — после обслуживания заявки с
преимуществом линия возвращается на дообслуживание прежней
заявки.
Мы рассмотрели ситуации, для которых характерна
очередь заявок. На практике встречаются также случаи
очереди линий, т. е. поступившая в систему заявка
находит свободные линии -и занимает одну из них в
соответствии со специальными правилами, аналогичными
правилам для очереди заявок (,в порядке освобождения линяй,
в случайном порядке, с учетом равномерной загрузки
линий и т. д.).
Для удобства изучения реальный (процесс
обслуживания можно представить в виде последовательности
различных фаз. Это равносильно расчленению системы
па автономные агрегаты таким образом, чтобы
последующий агрегат мог приступить к обслуживанию заявки
лишь после того, как работа предыдущего агрегата с
данлой заявкой полностью закончена. Примером
трехфазной системы может служить обслуживание
покупателей в -маг'азине: первая фаза — выбор товаров и
выписывание чека, вторая фаза — оплата чека в кассе, третья
фаза — получение покупки в отделе контроля и выдачи
покупок. Более сложным примером многофазного
обслуживания может быть технологический «процесс,
связанный с обработкой деталей на многостаночной линии.
Деталь не может быть принята к обслуживанию /г+1-м
агрегатом (станком), пока не будет закончена вся
работа, связанная с я-м агрегатом (станком).
Итак, мы рассмотрели наиболее распространенные
виды процессов функционирования систем массового
обслуживания. В заключение необходимо сказать иесколь-
63

ко слов об искомых величинах при решении задач,
связанных с системами массового обслуживания.
Для систем с отказами наиболее характерными
показателями качества обслуживания является средняя доля
отказов R(t0, t) за время (/0, to + t). Эта величина
определяется следующим образам. Количество заявок,
поступивших на обслуживание за время (/0, to + t), обозначим
Л/(/о, t), а среднее количество отказов за то же время —
mAoJ), тогда
Кроме средней доли отказов, для характеристики
качества обслуживания часто пользуются вероятностью
того, что за время (to,to + t) не будет ни одного отказа
P(tt t)
)
Для систем с ожиданием характеристиками качества
обслуживания обычно являются среднее время
ожидания, среднее количество заявок в очереди и т. д.
Качество обслуживания системами смешанного типа
характеризуется как средней долей отказов (или
вероятностью того, что за данный промежуток времени все
заявки будут обслужены), так и средним временем
ожидания (средним количеством зая.вок в очереди).
14. Моделирование потоков заявок
Рассмотрим способы моделирования на ЭВМ потоков
заявок, поступающих ,в систему массового обслуживания.
Сначала остановимся на достаточно простом и вместе с
тем наиболее распространенном случае, когща в систему
поступает ординарный стационарный поток однородных
событий с ограниченным последействием (поток типа
Пальма).
Любой поток типа Пальма может быть задан
функцией плотности f(z) случайных интервалов zy, / = 2, 3,...
между последовательными моментами /j = //_i+27-
поступления заявок. Для моделирования его на ЭВМ
достаточно построить необходимое число реализаций потока, т. е.
таких неслучайных последовательностей /ь t2y..., момен-
64

тов поступления заявок в систему, интервалы между
которыми Zj~tj—/,_i являлись бы возможными
значениями случайных величин Zj, описываемых функцией
плотности Hz).
* *
Процедура построения последовательности tu t2, ...
состоит в следующем. Сначала формируется t\ = Z\.
Функция плотности f\{z\) может быть определена через /(г)
по формуле Пальма D.4). Тогда, исходя из наличия в
ЭВМ электронного или алгоритмического датчика
случайных чисел Хг с равномерным распределением в
интервале @,1), приступаем к формированию z\. Для этого
одним из способов, рассмотренных в главе П,
преобразуем случайное число Х{ в случайное число Z\, имеющее
функцию плотности f\{z\). Получив гь полагаем
/Г= 7Х • D.8)
Далее переходим к формированию /2=^1+22.
Очередное случайное число х% аналогичным образом
преобразуем в случайное число г2, имеющее функцию плотности
f(z). Получив z2y полагаем
D.9)
В дальнейшем процедура получения любого tj
совпадает с процедурой формирования t2, т. е. очередное
случайное число преобразуется в случайное число Zj,
имеющее функцию плотности f(z) и
tl-tU + Ъ • D.10)
Рассмотрим некоторые примеры, часто
встречающиеся при решении практических задач методом
статистического моделирования.
Пример 1. Простейший (пауссоновский) поток.
Как отмечалось выше, простейший (пауссоновский)
поток является стационарным ординарным потоком
однородных событий без последействия, т. е. одним из
возможных частных случаев потоков типа Пальма. Функция
плотности /(г) для простейшего потока имеет вид
показательного распределения D.6):
5. Заказ 3012. 65

где X — интенсивность потока, определяющая среднее
значение числа заявок, поступающих в единицу времени.
Чтобы определить функцию плотности f{ (z\) для
первого интервала, воспользуемся формулой Пальма.
/iBri)=X[l—Я S e-^dz] • D.11)
о
После несложных вычислений получаем:
fi(zi)=berb± . D.12)
Отсюда следует, что функция плотности .первого
интервала для простейшего потока имеет тот же вид, что и
f(z). Этим свойством в общем случае, не обладают
другие потоки типа Пальма.
Таким образом, для формирования реализаций
простейшего потока необходимо иметь последовательность
случайных чисел z\, z2,..., имеющих показательное
распределение D.6) с параметро.м К. Методику получения
такой последовательности мы рассматривали в главе II.
В соответствии с B.15) случайные числа Zj могут быть
получаны по формуле
— 1
Zj=—-In (*<).
где Xi — случайные числа с равномерным распределением
в интервале @,1).
*
Последовательность моментов поступления заявок tj
будет следующей:
D.13)
Пример 2. Поток с равномерным распределением
интервалов Zj.
Функция плотности рассматриваемого потока имеет
вид равномерного распределения:
— для 0<z<b
D 14)
0 для г^О; z^b K }
66

Можно показать, что среднее значение
(математическое ожидание) случайной величины z равно-к-. Поэтому
среднее число заявок, поступивших в единицу .времени
(интенсивность потока):
*=f- .D.15)
Определим функцию плотности ]\{zx) для первого
интервала. По формуле Пальма аналогично формуле D.11):
D.16)
о
Заметим, что среднее значение длительности первого
интервала может быть получено как математическое
ожидание случайной величины, имеющей функцию
плотности D.16):
ЩгД~. D.17)
з
Перейдем к формированию первого интервала z{. Для
этого, имея случайные числа хг- с равномерным законом
распределения в интервале @,1), необходимо получить
случайное число zu соответствующее функции плотности
D.16). Преобразуем соотношение D.16) следующим
образом. Подставим в него вместо величины b ее значе-
2
ние Ь=—из равенства D.15). Тогда можно записать
функцию плотности
A=-f- D.18)
Легко усмотреть, что f\(z\) в таком виде совпадает с
функцией плотности B.16). Поэтому для определения zx
можно воспользоваться формулой B.20):
7у7,). D.19)
67

Формирование Zj\ /=2, 3,... сводится к получению
случайных чисел, имеющих равномерное распределение
в интервале (о, Ь) соответствующих функций плотности
D.14). Для этого достаточно случайные числа Х\ с
равномерным распределением в интервале @,1) привести к
интервалу (о, 1з):
7г = Ьх{. D.20)
Последовательность моментов поступления заявок
рассматриваемого потока определяется системой
уравнений D.13).
Аналогично могут быть описаны процедуры
формирования других потоков однородных событий с
ограниченным последействием. Соответствующие соотношения
приводятся в работах [3], [4] и др.
Необходимо отметить, что потоки, рассмотренные в
примерах 1 и 2, отличаются благоприятной особенностью:
интегралы в формуле Пальма и формулах
преобразования случайных чисел берутся в конечном виде. В общем
случае эти интегралы могут оказаться неберущимися.
Кроме того, функции плотности f(x) иногда задаются
таблично по результатам обработки статистического
материала. На практике в такого рода ситуациях
пользуются приближенными методами. Интеграл в формуле
Пальма вычисляется обычно для заданного набора Z\
численными методами. Это не оказывается существенно
на объеме вычислений, так как формула Пальма
применяется только один раз для данного потока.
Преобразование случайных чисел выполняется, как правило,
методом кусочной аппроксимации функции плотности в
соответствии с соотношениями B.23), B.24) и B.25).
15. Моделирование одноканальной системы
массового обслуживания
Изучение методики моделирования систем массового
обслуживания начнем со случая одноканальной системы.
Итак, в систему поступает ординарный поток заявок
с заданным законом распределения. Предельное время
ожидания заявки начала обслуживания т—случайная
величина с законом распределения ф(т)» Время занято-
68

сти канала (длительность обслуживания) задается
законом распределения /(т*).
Целью моделирования является получение
характеристик качества обслуживания: среднего времени
ожидания в очереди, доли обслуженных заявок, доли заявок,
получивших отказ .и т. д.
Процесс функционирования системы рассматривается
за период времени [О, Г], т. е. заявки, для которых момент
появления tj>T, в систему не попадают и не
обслуживаются. Заявки, для которых время окончания
обслуживания больше Т, считаются (Получившими отказ.
Для построения алгоритма, моделирующего процесс
функционирования одноканальной .системы массового
обслуживания, нам понадобятся следующие операторы:
Ф\ — формирование очередного момента tj
поступления заявки в систему;
Р2 — проверка условия lj<T принадлежности
момента поступления очередной заявки tj интервалу
исследования системы [0, 71;
* *
Р3 — проверка условия tj<tj-\, где tj-\ — момент
окончания обслуживания предыдущей заявки;
Ф4 — формирование предельной длительности
ожидания заявки до начала обслуживания т в соответствии с
законом распределения ф (т);
As — вычисление момента t/ = tj + x (в момент t/
заявка покидает систему, если она не будет принята к
обслуживанию);
Ре — проверка условия t/<t^\ (заявка покидает
систему ранее чем освободится канал);
F7 — выбор в качестве момента начала обслуживания
/-й заявки момента окончания обслуживания /—1-й
заявки //Н)=/Д;
Fg — выбор в качестве момента начала обслуживания
/-й заявки момента ее поступления в систему ?jH)=/;-;
Ф9 — формирование длительности обслуживания т*
заявки (времени занятости канала) в соответствии с
законом распределения /(т*);
Аю — вычисление момента окончания обслуживания
/-й заявки (момента освобождения канала) tj = tjH) +т*;
Ри—проверка условия /* =^Г;
К\2 — подсчет количества обслуженных заявок т;
69

Ai3 — вычисление длительности пребывания в
очереди для /-й заявки tj—/j(H);
Кн — подсчет количества заявок п, получивших отказ;
К15 — подсчет количества (реализаций N;
Pie — проверка уславия N<N*> где #*—количество
реализаций, необходимое для обеспечения заданной
точности;
Fi7 —-переход к очередной реализации;
Aia —обработка результатов моделирования;
Ям — выдача результатов.
В операторной форме моделирующий алгоритм для
процесса функционирования одноканальяой системы
массового обслуживания имеет следующий вид:
13'14'17Ф1 Р2|15 Рзц Ф4
А5 Ре" F? SF87'8 Ф9 D.21)
AioPn;i4Ki2 AiJ 6>liKi4
2Ki5 Pie;i8 тЪ^АпЯм.
Исходными данными для моделирования являются
граница интервала .времени Т, законы распределения
ср(т) и f(x*)9 количество реализаций, обеспечивающих
заданную точность расчета N*, а также состояние
системы при*/ = 0, например /j_i = 0; /п = 0; п = 0\ N = 0.
Для наглядности изображения алгоритма приводим
его блок-схему (см. рис. 8).
Остановимся кратко на работе алгоритма и
отдельных его операторов. Оператор Ф1 формирует случайные
числа (интервалы времени между заявками) в
соответствии с законом распределения потока заявок.
Оператор Рг проверяет, принадлежит ли данная
заявка интервалу «времени [О, Т]. Если нет (t^T)
(реализация закончилась), то по стрелке с индексом 0 управление
передается оператору Kis Для подсчета количества
реализаций. Бели условие, проверяемое оператором Р2,
выполнено (tj<T) —заявка принадлежит рассматриваемой
реализации; тогда управление передается оператору Рз.
Этот оператор посредством проверки неравенства
ti<tj-i определяет, свободен или занят канал
обслуживания. Если tj^tj-\ (канал свободен), то. по стрелке с
70

Счетчии числа об'слд-
же/шь/х зплдок т
t{H)- t
7 7
Счегчин числа ошазодп
Переход и очередной
реализации
Обработка результо -
то8 моделирования
Выдача результатод
Рис. 8.

индексом 0 управление передается оператору F8, который
формирует i(p —момент начала обслуживания /-й
заявки (в данном случае tjH) =tj).
Если канал занят обслуживанием предыдущей заяв-
*
ки, tj<tj-\, по стрелке с индексом 1 управление
передается оператору Ф4, который формирует tj — время ожидания
в очереди — и передает управление оператору А5.
Оператор As вычисляет tj' = tj+Tj— момент, когда заявка
получит отказ, если не будет принята к обслуживанию, и
передает управление оператору Рб, который проверяет,
освободится ли канал к моменту t/. Если //</j_b т. е.
канал занят и заявка не может быть обслужена, тогда по
стрелке с индексом 1 управление передается оператору
Кн, который фиксирует число п заявок, получивших
отказ, и передает управление оператору Ф\ для
формирования tj+\. Если к моменту tj' = tj+tj, канал будет
свободен (tj'^tj-i), управление по стрелке с индексом 0
передается оператору F7. Оператор FV формирует момент
начала обслуживания /-й заявки (в данном случае это
будет момент освобождения канала ^(н) =^-0 и
передает управление оператору Ф9, формирующему Tj* —
длительность обслуживания /-й заявки. Оператор Аю
вычисляет tj* — момент окончания обслуживания (/j* = ^(H) +
+ Ti*) — и передает управление оператору Рц для
проверки неравенства tj*^T. Если это неравенство
выполнено, т. е. момент окончания обслуживания принадлежит
данной реализации, управление по стрелке с индексом 1
передается оператору Кп для подсчета количества т —
обслуженных заявок, а затем оператору А^ для
вычисления tj(H)—tj—'времени ожидания заявки до начала
обслуживания. Если неравенство tj*^T не выполнено, т.е.
момент окончания обслуживания tj* находится вне
данной реализации, управление передается оператору Кн,
учитывающему случаи отказа. Группа операторов
Ki5... Я19 играет в алгоритме вспомогательную роль и не
связана непосредственно с моделированием процесса.
Оператор Kis подсчитывает количество реализаций N.
Оператор Pie сравнивает N с Af*, где jV* — требуемое
количество реализаций. Если N недостаточно велико
(N<N*), управление по стрелке с индексом 1
передается оператору Fi7, осуществляющему переход к очередной
12

реализации. Если N = N* (реализаций достаточно)—по
стрелке с индексом 0 управление передается оператору
Ais для обработки результатов моделирования
(обработка данных, накопленных операторами: Ki2 —
количество обслуженных заявок, Ки — количество отказов,
Ai3 — длительность -пребывания заявки в очереди) и
затем оператору Яю для выдачи этих результатов на
печать.
Мы подробно остановились на работе алгоритма,
моделирующего одноканальную систему MaocoiBoro
обслуживания, в связи с тем, что одноканальная система,
являясь простейшей из всех видов систем массового
обслуживания, тем не менее с методической точки з;рения содер-
Ж'ит многие элементы, свойственные системам более
общего характера.
16. Моделирование многоканальной системы
Рассмотрим теперь структуру моделирующего
алгоритма для случая многоканальной системы.
Пусть система массового обслуживания имеет п
каналов, вполне идентичных (с точки зрения длительности
обслуживания) и равноправных (в смысле возможности
использования для обслуживания данной заявки любого
свободного канала).
Будем ло-прежнему считать, что в систему
поступает ординарный поток заявок с заданным законом
распределения. Время возможного ожидания до начала
обслуживания т — случайная величина с законом
распределения ф(т). Длительность обслуживания т* имеет закон
распределения /(т*). Процесс функционирования
системы рассматривается в интервале времени (о, Т).
Заявки принимаются к обслуживанию в порядке очереди (в
порядке поступления в систему). Если для
обслуживания данной заявки имеется несколько свободных
каналов, то каналы привлекаются к обслуживанию в порядке
очереди (в первую очередь привлекается тот канал,
который ранее других освободился от обслуживания).
Описанная здесь система близка по характеру к одно-
канальной системе. Поэтому многие операторы схемы
D. 21) могут быть использованы для построения
моделирующего алгоритма в нашем случае. Сохраним также
73

все обозначения, использованные для построения схемы
D.21).
Рассмотрим следующие операторы:
Аю — вычисление момента окончания обслуживания
/-й заявки (момента освобождения k-то канала) t*K =
= 'j^ +т*, где t$ — момент начала обслуживания 1-й
заявки /е-м каналом;
А20 — запись полученного VlK в регистр;
A2i — выбор min tiK из регистра — наименьшего из
каналов, свободных в данный момент времени;
Р3 — проверка условия /j<min t*K (в момент
поступления заявки все каналы заняты);
Рб —.проверка условия l'j<min tfK (заявка покидает
систему ранее чем освободится хотя бы один из
каналов);
F7— выбор в качестве момента начала
обслуживания Кн) величины min i*K ;
Рп — проверка условия tt*^.T (период пребывания
t-н заявки в системе целиком принадлежит интервалу
(о; Т) моделирования процесса).
Запишем операторную схему моделирующего
алгоритма для многоканальной системы, пользуясь также
операторами схемы D.21):
13,14,17ф1 Р121 2ф3|8
2 jib
Ф4 А5 Р,^14 F? D.22)
3F8 7'8Ф9 А?о 20Piui4 Ki2
Ajs 641Kl4 2Kl5 P16Ц8 FJt
16A18 Я» 10A2J 2A2Si.
Ha pnc. 9 представлена блок-схема рассматриваемого
алгоритма.
Моделирующий алгоритм D.22) работает аналогично
алгоритму D.21). Поэтому мы не будет
останавливаться подробно на его описании. Отметим только некоторые
особенности.
74

Рормиро-
6ание z
i—L
Формирование tj
Выбор минимального ten из регистра
W
20
Формирование с *
_L
+ г
Запись полученного tpH 8 регистр
П
JL
Счетчик числа обслуженных за л вон т
±
tj~i
Счет чин числа отназод /7
/V+/
17
Переход н очередной реализации
18
Обработна результатов моделирования].
19
I
Выдача результатов
Рис. 9.

Когда оператором Аю определены момент tfK
окончания обслуживания 1-й заявки ?-м каналом, имеется, с
сдной стороны, момент окончания обслуживания 1-й
заявки, а с другой — момент освобождения &-го канала.
Далее t*K , как момент окончания обслуживания /-и
заявки, сравнивается с T (оператор Рц). Если условие,
(проверяемое оператором Рц, не выполнено — заявка не
считается обслуженной и управление передается
оператору Кн Для подсчета числа отказов. Когда же условие,
проверяемое оператором Рц, выполнено, управление
передается оператору Ki2 для подсчета числа обслуженных
заявок. Кроме того, t*K рассматривается в этом случае
как момент освобождения &-го канала. Может случиться,
что в данный момент времени имеются и другие
свободные каналы. Моменты освобождения их помещены -в
специальный регистр (ячейку памяти ЭВМ) в порядке
возрастания. Оператор Аго записывает в этот регистр
очередное значение t*K , полученное оператором Аю. В
отличие от схемы D. 21) от оператора Рг по стрелке с
индексом 1 управление передается оператору А2ь который
выбирает из регистра наименьшее t*K и передает
управление Рз. Таким образом, наличие свободного канала
выясняется сравнением tj — момента поступления
очередной заявки —с riling* —минимальным временем
освобождения одного из каналов. В случае
отрицательного решения вопроса о наличии свободного канала, в
момент поступления заявки, величина mint*K
используется (оператор Ре) для определения возможности
обслужить заявку. Если данная заявка может быть
обслужена, то в качестве момента начала обслуживания /^н)
опять же выбирается min /,* (оператор F7). В
остальном работа алгоритма D.22) не отличается от работы
алгоритма D. 21).
В заключение заметим, что структура
моделирующего алгоритма существенно зависит от порядка выбора
заявок из очереди для обслуживания (дисциплина
очереди заявок) и порядка привлечения к обслуживанию
свободных каналов (дисциплина очереди каналов).
Некоторые наиболее распространенные варианты
дисциплины очереди заявок и каналов мы рассмотрим в
следующем «параграфе.
76

17. Модификаций систем массового обслуживания
Системы массового обслуживания как
математические схемы для формального описания
функционирования реальных систем крайне дифференцированы.
Поэтому случаи, встречающиеся на практике, далеко
не исчерпываются рассмотренными выше простейшими
одноканальной и многоканальной системами. Виды этих
«систем столь многочисленны, что нет возможности даже
.кратко на всех остановиться. Теперь мы рассмотрим
важнейшие виды систем массового обслуживания более
общего характера, чем описанные выше; приведем
соответствующие примеры.
Начнем с возможных вариантов дисциплины очереди
заявок. В рассмотренных выше одноканальной и
многоканальной системах массового обслуживания заявки,
поступающие в систему, обслуживаются в порядке
поступления. Иногда для некоторых важных в
практическом отношении задач такой порядок обслуживания
оказывается невыгодным. Нередки случаи, когда заявки,
имеющие право (в порядке очереди) на немедленное
обслуживание, могут находиться в системе довольно
долго (т — велико) и обслуживаться после других заявок,
имеющих малое т и поэтому быстро покидающих
систему. Очевидно, что здесь мы сталкиваемся с упомянутым
выше случаем целесообразности изменения порядка
обслуживания, а именно: предпочтительного обслуживания
заявок, ранее других покидающих систему.
Рассмотрим несложную модификацию, позволяющую
использовать алгоритм D.21) для моделирования
систем такого типа.
В момент освобождения линии t*j-\ в очереди может
оказаться несколько заявок. Необходимо предусмотреть
специальный блок операторов (подалгоритм),
занимающихся фиксацией всех // для заявок, находящихся, в
очереди, и отысканием заявки, которая ранее других
покидает систему.
Для модификации алгоритма D.21) нам требуются
следующие операторы:
Pi—проверка условия k>\ (где k — количество
заявок в очереди);
А2 — определение /'mm — наименьшего значения /'
77

5
!
Счетчик числа
назов л
у
М-1
от-
(момент окончания
ожидания) для заявок,
находящихся в очереди;
Ръ — проверка условия
//min<^-i (совпадает с
оператором Рб в
моделирующем^ алгоритме D.21);
К4 — счетчик количества
п — заявок, получивших
отказ вследствие
недостаточного времени ожидания
(совпадает с Км);
Ks — счетчик количества
заявок в очереди
(вычисляет k— 1);
Рб — проверка условия
/г>0 (где k — количество
заявок в опереди).
Запишем операторную
схему этого подалгоритма:
Кь
D.23)
Блок-схема его
приведена на рис. 10.
Порядок работы
подалгоритма следующий.
Если в очереди имеется
k заявок (проверка этого
факта осуществляется оператором Pi (?>1), надо
сначала выбрать заявку с наименьшим V — моментом
окончания ожидания. Этим занимается оператор А2.
Определив fmin, A2 передает управление Рз, который
проверяет, может ли данная заявка быть обслуженной
(t'min^t*^). Если ^min^f/.p управление передается
оператору, формирующему #н> — момент начала
обслуживания (/(H) = /J_1). Если заявка обслужена не будет,
то через оператор K-ь подсчитывающий п — количество
отказов, управление передается оператору As для вычис-
78

ления количества заявок в очереди (k—1) и затем
оператору Рб для проверки неравенства /г>0. Если k>0
(существует очередь заявок), управление передается снова
оператору Pi и затем А2 для перегруппировки заявок в
очереди (определения нового /'mm). Если неравенство
&>0 не выполнено (нет очереди заявок), управление
передается оператору, формирующему /j+j — момент
поступления новой заявки.
Эта модификация моделирующего алгоритма
позволяет кроме обычных результатов получить
дополнительные характеристики моделируемого процесса, такие,
как среднюю длину очереди и среднее время ожидания.
Это нетрудно сделать, выводя на печать величину k для
каждого момента изменения его значения (поступления
новой заявки, взятие очередной заявки на обслуживание
или получения отказа), а также моменты t'j окончания
ожидания для всех заявок >в очереди.
В качестве примера других возможных вариантов
дисциплины очереди рассмотрим обслуживание в
случайном порядке -в соответствии с заданными
вероятностями.
Пусть время пребывания заявки в системе
ограничено жак в алгоритме D.21), а заявки принимаются к
обслуживанию в соответствии с вероятностями qa,
задаваемыми заранее или вычисляемыми в процессе
моделирования.
Приводим необходимые операторы:
Рп —проверка условия fj<t*»x\
Kn+i — счетчик числа отказов га;
Ап+2 — запоминание (запись ,в регистр) заявки,
удовлетворяющей условию; проверяемому
оператором Р„;
К?м-з — счетчик количества у — заявок в регистре
(операция Y+Oi
Кп+4 — счетчик количества k — заявок в очереди (k—1);
Рп+5 — проверка наличия заявок в очереди (&0)
Рп+6 — проверка наличия заявок в регистре
Ап+7 — выбор заявки в соответствии с заданными
вероятностями;
А?1+8 — формирование #н) — момента начала
обслуживания для выбранной заявки (tw = tj_{).
79

Операторную схему подалгаритма можно записать
так:
я+1
К
и
n+4
Pn
Of
пч
Счетчин
п+г
Зопис
л>3
Ч"н у
У1
число о/пнозов п
ь t j в регистр
\
\
Вы fop заядни по
жребию
п+8 \
J-'
Кп+3
D.24)
+6
Рис. 11.
На рис. 11 приводится
блох-схема
подалгоритма. Рассмотрим кратко
его работу.
Оператор Рп с
помощью цикла,
контролируемого операторами
Кп+4 и Рп+6, отбирает
заявки, которые могут быть
приняты к
обслуживанию, т. е. заявки, для
которых t'j>tj-u моменты
времени tj для этих
заявок заносятся в регистр
Ап+2 (из них
впоследствии и будет выбрана
заявка в соответствии со
жребием). Когда все
заявки в очереди
просмотрены, от оператора Рп+5
по стрелке с индексом О
управление передается
оператору Яп+б, который
определяет наличие
заявок в регистре и, если
они там имеются (y>0)>
оператор An+z выбирает
одну из них в
соответствии с заданными
вероятностями. Затем для этой
заявки оператор An+s
формирует момент
начала обслуживания. Если
же заявок в регистре нет,
а это значит, что ни одна
заявка из очереди не в со-
80

стоянии дождаться момента освобождения канала и,
следовательно, все они получили отказ (что было
зафиксировано оператором Rn+i), тогда управление от
оператора Рп+6 по стрелке с индексом 0 передается оператору,
формирующему новую заявку.
Приведенная блок-схема подалгоритма может быть
попользовала для моделирования простейшего случая
обслуживания с преимуществом, а именно для выбора
заявки в (соответствии с заданным коэффициентом
преимущества.
Действительно, если заменить оператор Ап+7 (выбор
заявки по жребию) оператором — .выбор заявки в
соответствии с коэффициентом преимущества, мы и получим
искомый алгоритм.
Переходя к рассмотрению дисциплины очереди
обслуживающих каналов многоканальной системы массового
обслуживания, остановимся ий алгоритме,
моделирующем выбор канала в соответствии с заданным правилом.
Для этого нам потребуются следующие операторы:
Рг —проверка условия i>n, где i — текущее, an —
максимальное значение номера обслуживающего
канала;
Кг+1 — счетчик числа каналов (реализуетоперацики+1);
Рг+2 — проверка условия ^</}св)где tj — момент
поступления очередной заявки, a tK?0)— момент
освобождения f-ro канала;
Fr+3 — запоминание (занесение в регистр) ^с),для
которых не справедливо условие, проверяемое
оператором Рг+2*,
Кг+4— счетчик количества свободных каналов в регистре
(реализует операцию я*+1, где /г*— количество
свободных каналов);
Рг+5 — проверка условия п*>0 (т. е. имеются ли
свободные каналы в регистре);
Fr+6 — реализация заданного правила выбора канала для
обслуживания.
Операторная схема подалгоритма имеет вид:
D.25)
Kr+4 rP
Блок-схема приведена на рис. 12.
Подалгоритм работает следующим образом.
6 Заказ 3012.

X
г+3
Запись t
г + Ч
(сб)
1
^
/Г-
>
г )
,0
8 регистр
2 -*
г+6
Реализация правило
бы бора нонало
и
Рис. 12.
Для каждой
поступившей в момент tj заявки
перебираются все
обслуживающие каналы и
выявляются свободные.
Цикл перебора
обеспечивается операторами Рг и
Кгч 1. Занятость капала
оценивается посредством
Рг-ь2. Незанятые каналы
фиксируются в регистре
(оператор Fr+3), а
количество их — в счетчике
Кг-н. Отобрав все свобод-,
ные к моменту tj каналы,
переходим через
оператор Рг+5 (л*>0) к опера-
тору Fr+6, реализующему
правило выбора канала.
От оператора Fr+5
управление передается
оператору, 'формирующему //н>
— момент начала
обслуживания /-й заявки
выбранным каналом. Если
условие, проверяемое
оператором Рг+5, не
выполнено (л* = 0), это
значит, что мы, перебрав все
каналы, не обнаружили
(а потом и не записали в
регистр) ни одного
свободного к моменту tj канала; тогда мы передаем
управление оператору, формирующему новые значения tj и
отбираем свободные каналы для момента /j+i-
Рассмотренные здесь варианты дисциплины очереди
и порядка обслуживания для одноканальнои и
многоканальной систем массового обслуживания понадобятся
нам в дальнейшем для моделирования реальных систем
массового обслуживания.
Следующим важным обобщением, требующим
специального рассмотрения, являются системы массового
обслуживания с ненадежными элементами.
82

Во многих практических задачах системы массового
обслуживания рассматриваются как системы с
ненадежными элементами. В общем случае это означает, что
в некоторый момент времени/(сб) в работе
обслуживающего канала происходит отказ (сбой); канал может быть
отремонтирован за время т(р) и снова приступить к
обслуживанию в момент времени t = t{QC) +т(р).
Надежность канала задается функцией распределения
интервала безотказной работы P(t)> равной вероятности
того, что /<сб)>/:
P{t)=P{t^)>t)y D.26)
а тA)) (время ремонта)—как случайная величина с
заданным законом распределения.
Можно интерпретировать сбой и ремонт канала как
обслуживание фиктивной заявки с моментом начала
обслуживания Исб\ длительностью обслуживания ^ и
главенствующим коэффициентом преимущества. Весь
период ремонта канал можно считать фиктивно занятым.
Обычно принимаются два предположения относительно
судьбы заявки, находящейся на обслуживании в момент
сбоя; либо заявка получает отказ и уходит из системы,
либо ее обслуживание может быть продолжено после
устранения последствий сбоя (после ремонта).
Рассмотрим принципы моделирования систем с
ненадежными элементами на примере одноканалыной
системы D.21).
-Пусть мы имеем одноканальную систему. Дополним
ее описание следующими свойствами: в моменты времени
t/c6\ заданные законом распределения P(t), канал
выходит из строя; ремонт его продолжается случайное время
т(р>, заданное законом распределения f(x); заявка, при
обслуживании которой наступил сбой, получает отказ.
Воспользуемся следующими операторами:
Р8 — проверка условия а>0;
ФбЧ.! — формирование случайного числа t{c6) в
соответствии с функцией распределения D.26);
Fs+2 — выбор /(сб) из памяти;
Ps+3 — проверка условия ^(сб)</(сб);
Fs+4 — запись /(сб) в память, формирование а+1;
F3+5 — формирование а = 0;
Ps+6 — проверка условия /<сб><*п;
Фй+7 — формирование случайного значения т<р) —
времени ремонта обслуживающего канала;
6* 83

As+8 — определение момента №— готовности капала к
обслуживанию после ремонта;
Р5+9 — проверка условия #г)<#ж) (возможности
обслужить заявку после ремонта).
Операторная схема подалгоритма учета надежности
может быть представлена >в виде:
Блок-схема подалгоритма представлена па рис. 13.
Подалгоритм работает следующим образом.
Оператор Ps проверяет условие а>0. Если а>0, это
свидетельствует, что значение t{c6) уже сформировано и
записано в память; управление передается оператору Ф«+2
для извлечения /(сб> из памяти. В противном случае
(а = 0) оператор Fs+i формирует значение №б\ Оператор
Ps+з проверяет неравенство /<сб)<^сб). Если это
неравенство не выполнено, т. е. сбой произошел после окончания
обслуживания, то заявка может быть обслужена. В этом
случае значение tt{c6) следует записать в память,
сформировав признак а=1 (свидетельствующий о наличии ?(сб)
в памяти). Если неравенство /(сб)<#св> выполнено,
вопрос о возможности обслужить заявку будет
рассматриваться особо. Значение /(сб)уже не понадобится нам прк
рассмотрении последующих заявок, поэтому управление
по стрелке с индексом 1 передается оператору Fs+5 Для
формирования а = 0. Если неравенство /<сб)<#н),
проверяемое оператором Ps+6, не выполнено, т. е. #сб)>Лн)
(сбой произошел после начала обслуживания)—заявка
получает отказ (управление по стрелке с индексом 0
передается оператору, фиксирующему количество отказов).
Если условие, проверяемое операторам Ps+6,
выполнено— /(сб)</(ц) (сбой произошел до начала
обслуживания), вопрос о возможности обслужить заявку нельзя
.решить, не зная т(р)—времени ремонта. Оператор Ф<?+7
формирует т(р) и затем оператор As+8 вычисляет /(г) = /(сб) +
+%М (момент окончания ремонта). Тетерь судьба заявки
может быть выяснена. Если условие, проверяемое
оператором Ps+9, выполнено — /<г></<>«> (ремонт окончится
раньше, чем истечет время ожидания), управление по
84

стрелке с индексом 1
передается группе
операторов,
занимающихся обслуживанием
заявки; если условие,
проверяемое
оператором Ps+9, не
выполнено (заявка не
дождется окончания
ремонта) — по стрелке с
индексом 0 счетчлку
количество отказов п.
Аналогично могут
быть построены
моделирующие алгоритмы
для более сложных
случаев учета
надежности каналов системы
массового
обслуживания (например,
зависимость следствий
отказа от момента его
наступления и др.)-
Дальнейшие
обобщения обычно связаны
с учетом
неоднородности заявок. Каждая
заявка
.характеризуется некоторым набором
параметров,
описывающих ее
индивидуальные свойства. Это
обстоятельство требует
особых методов
описания и моделирования
потока заявок и
процесса
функционирования самой системы
массового облужива-
ния.
Запись t(c^ в память
Формирование а=/
Рис. 13.
Пусть /-я заявка описывается случайным моментом
t
поступления в систему tj и параметрами
о&2>,
85

принимающими случайные значения. Математической
схемой для описания потока неоднородных заявок
может служить схема потока случайных векторов. В общем
виде эта схема сложна и громоздка с точки зрения
реализации на ЭВМ. При решении практических задач
стремятся поток моментов поступления заявок представить
как поток однородных событий, а случайные параметры
описать соответствующими условными (или
безусловными, когща параметры и поток моментов независимы)
распределениями вероятностей. При этом мы приходим
к необходимости формировать па ЭВМ реализации
многомерных случайных векторов (см. главу II). Объем
вычислений уменьшается, если составляющие
случайного вектора независимы (например, поток с
ограниченным последействием).
Неоднородность заявок оказывает влияние и на
процесс «их обслуживания. Наиболее распространены
случаи, «когда характеристики обслуживания, такие, ,как
время ожидания обслуживания т, время занятости канала
(длительность обслуживания) т* и т. д., являются
функциями параметров заявки аь «2, . . . , а&. Это
обстоятельство учитывается при моделировании: к операторам
моделирующего алгоритма добавляются соотношения
для расчета % и т*.
Наконец, на практике встречаются системы,
обладающие специальными комплексами управления, которые
по заданным параметрам заявки и известным
характеристикам состояний системы обслуживания
вырабатывают команды, определяющие порядок обслуживания и
значения характеристик обслуживания. Моделирование
таких систем на ЭВМ проводится с использованием
рассмотренных выше методов, но с дополнениями к
моделирующему алгоритму, имитирующими работу комплекса
управления. Обычно имитация сводится к
воспроизведению на ЭВМ некоторых формально описанных
алгоритмов, хотя бы приближенно представляющих работу
комплекса управления.

Глава V
МОДЕЛИРОВАНИЕ АВТОМОБИЛЬНОГО
ДВИЖЕНИЯ
18. Вводные замечания. Формальное описание дороги
Опыт построения моделирующих алгоритмов
для систем массового обслуживания оказывается
полезным при моделировании многих типов сложных систем,
встречающихся на практике. Наиболее близки по
структуре и приемам моделирования к процессам массового
обслуживания различные виды производственных
процессов. Для формального описания производственный
процесс представляется в виде последовательности
производственных операций, главнейшими из которых являются
операции обработки деталей (полуфабрикатов),
сборки изделий, операции управления технологией и т. д.
Каждая операция может быть описана как некоторый
процесс массового обслуживания. Например, операция
обработки детали есть «обслуживание» ее
соответствующим станком. Таким образом, производственный процесс
может быть представлен в виде многофазной системы
массового обслуживания. Методика моделирования
производственных процессов широко описана в литературе
по статистическому моделированию (см. например,
работы [2], [3], [4] и др.), поэтому мы на ней
останавливаться не будем.
Мы изложим основные черты методшш
моделирования процесса автомобильного движения, существенно
отличающегося от процессов массового обслуживания.
Естественно, что некоторые понятия, лежащие в
основе формального описания процессов функционирования
систем массового обслуживания, и методические приемы
87

построения моделирующих алгоритмов для них могут
быть перенесены на автомобильное движение. Например,
автомобили, появляющиеся на дороге, можно
рассматривать как заявки; дорогу — как совокупность средств,
обеспечивающих обслуживание заявок. Однако, строго
говоря, автомобильное движение не может быть описано
в виде традиционной системы массового обслуживания.
Попытки попользовать методы теории массового
обслуживания для решения задач, связанных с исследованием
автомобильного движения, оказываются, как правило,
недостаточно перспективными. Основная причина
этого — особенности автомобильного движения как
процесса (зависимости характеристик обслуживания от
параметров заявок — типа автомобиля, его скорости и т. д. —
и от интенсивности движения на дороге — канал с
насыщением). Это обстоятельство приводит к необходимости
разработки специальной методики моделирования
автомобильного движения. Тем не менее опытом
моделирования процессов массового обслуживания мы будем
широко пользоваться в дальнейшем.
Сущность моделирования автомобильного движения
в предлагаемом ниже виде сводится главным образом к
определению состояний системы (координат всех
автомобилей и их скоростей) в некоторые моменты времени
A
Скорость автомобиля в данный момент времени,
определяемая балансом тяговой силы двигателя и
сопротивлением дороги, зависит от типа и состояния покрытия
дороги, от величины ,и направления уклонов продольного
профиля дороги и кривизны их в шлане. Кроме того,
скорость автомобиля может ограничиваться дорожными
знаками и правилами уличного движения (снижение
скорости на поворотах, перекрестках и т. д.). Поэтому при
формальном описании дороги необходима информация,
достаточная для вычисления скорости автомобиля в
любой момент времени. Обычно оказываются нужными
следующие параметры дороги: v — ширина покрытия дороги
(количество полос); / — коэффициент сопротивления
качения для покрытия данного типа; Д/ — поправка к / за
счет погодных условий (дождь, гололед ,и т. д.); i — уклон
продольного профиля дороги; \i — признак ограничения
скорости (за счет дорожного знака или поворота,
перекрестка и т. д.).

Для решения рассматриваемых ниже задач
достаточно следующее приближенное описание дороги.
Исследуемый участок дороги между точками А ,и В разбиваем на
интервалы, граничными точками которых являются так
называемые особые точки. К особым точкам относим
каждую точку, ib которой происходит изменение хотя бы
одного параметра дороги. Например, особыми точками
оказываются перевалы, повороты, границы зон действия
дорожных знаков, светофоры, точки изменения типа и
состояния покрытия дороги, а также точки изменения
уклона и
На реальных дорогах уклон профиля I изменяется
плавно. Точная запись величины i, как функции
расстояния, /привела бы к весьма громоздким соотношениям,
неудобным для ввода в ЭВМ. Поэтому зависимость i от
расстояния принимается приближенной и профиль
дороги представляется в виде кусочно-линейной функции с
постоянными .значениями i на отдельных участках.
Границы этих участков также являются особыми точками на
формализованной дороге.
Характеристики дороги задаются <в виде таблицы по
аргументу /—расстоянию данной особой точки от точки
А. Первой особой точкой всегда будет точка А. Величина
U — расстояние до следующей особой точки.
Таблица 3
Особые точки
(расстояние
в м)
i
f
А/
h
Л
АЛ
h
h
h
Л/а
1з
It
/з
Д/з
V-з
1п
In
fn
А/я
V-n
Для интервала длиной 1\ задаются (значения /j, f\, Af\
и jlii. Заметим, что в ряда случаев Д/ может быть равной
89

нулю, а значение [i отсутствовать. Далее задается
расстояние /2 до второй особой точки, а также
соответствующие значения параметров дороги /г, Ь, А/2, \i2 и т. д.
Естественно, что (некоторые параметры могут оставаться
одинаковыми на различных интервалах между особыми
точками.
При оценке влияния параметров дороги на скорость
данного автомобиля необходимо по его ^координате s
(расстоянию от точки А) войти в таблицу, определить
характеристики дороги и использовать их для вычисления
скорости автомобиля.
Такой метод представления данных о
характеристиках дороги достаточно точен для решения практических
задач ,и не требует для его реализации расхода памяти и
большого числа операций на ЭВМ.
Мы будем рассматривать только двухполосную
дорогу. Другие типы дорог описываются аналогично.
19. Описание потока автомобилей
Для решения задач, связанных с исследованием
процесса автомобильного движения, поток автомобилей
необходимо рассматривать как случайный поток с
заданными характеристиками.
Моменты появления автомобилей на дороге (в точке
А или в точке В) мы будем представлять как случайный
поток однородных событий в смысле теории массового
обслуживания. Исходя из характера движения
автомобилей можно считать этот поток стационарным или
нестационарным ординарным потоком с ограниченным
последействием. Поэтому интервалы времени между по-
следовательньши моментами появления автомобилей
оказываются независимыми случайными величинами с
соответствующими законами распределения.
Интенсивность движения на дороге будем
характеризовать интенсивностью потока автомобилей, движущихся
в одну сторону,— среднем числом X автомобилей,
проходящих в час через поперечное сечение дороги.
При подборе закона распределения интервалов
времени между моментами поступления на дорогу соседних
автомобилей нужно иметь в виду следующее. Обычные
90

законы распределения, применяемые ib теории массового
обслуживания (например, показательный, равномерный
и т. д.), допускают с некоторой вероятностью весьма
малые или близкие к нулю значения интервалов времени.
Правила уличного движения предусматривают, что
между соседними автомобилями должны выдерживаться
определенные интервалы (не менее /б), обеспечивающие
безопасность движения. Величина /б зависит от скорости
движения автомобиля и выражается некоторой функцией
скорости. Например, часто пользуются выражением
/б= 12 + 0,1 у1-85. E.1)
Поэтому на (практике моменты появления
автомобилей на дороге не могут быть сколь угодно близкими друг
к другу. Из практических соображений они должны
отличаться друг от друга не менее чем на некоторую
величину а. Тогда интервал .времени Zj между моментами
поступления автомобилей f,-_i и U могут быть выражены
суммой
Ъ = а + Ъг E.2)
где gj — случайная величина с заданным законам
распределения.
Рассмотрим некоторые примеры. Пусть поток
однородных событий характеризуется показательным
распределением интервалов z$
f(z)=X*e~X*Z E.3)
с интенсивностью Я*. Это значит, что мы имеем дело с
простейшим (|пуассоновским) потоком см. D.5), D.6).
Предположим теперь, что мы рассматриваем поток,
интервалы Zj которого имеют вид:
. E.4)
Тогда закон распределения для интервалов E.4)
будет выражаться функцией (плотности
для z ^а;
-%*{Г-а) E.5)
\*е для z>a. v ;
Интенсивность К потока E.5) можно определить по
формуле
91

Аналогично можно записать закон Распределения для
других случаев. Пусть поток характеризуется
равномерным распределением интервалов zy.
( 1
*/ \ "Г Для 0<г<6; E.7)
/(г)- ) Ь v '
( 0 для z^b\ г^О.
Легко видеть, что интенсивность такого потока в силу
равенства D.3)
о
* —
- т
Если теперь рассмотреть поток с интервалами
то закон распределения интервалов Zj будет иметь вид
равномерного распределения:
О для г^а\
— для a<z<b + a\ E.8)
О для z>b + a,
а интенсивность потока X определяется соотношением
Поток автомобилей «в общем случае состоит из
автомобилей п типов. Состав потока можно характеризовать
долей |3k (в %) от общей интенсивности потока X. Для
каждого типа автомобилей устанавливается
максимальная скорость Ущах как скорость свободного движения по
горизонтальной дороге без ограничений. Кроме того,
должен быть задан полный вес автомобиля с грузом. Вес
груза инопда считают случайной величиной с заданным
(для ,каждого типа автомобилей) законом распределения.
Фактическая скорость свободного движения
автомобиля с учетом профиля дороги, веса груза и тяги
двигателя в каждый момент времени определяется
дифференциальными уравнениями силового баланса. Для
определения скорости эти дифференциальные уравнения
92

должны быть решены. Однако 'решение
дифференциальных уравнений в процессе моделирования движения на
дороге для многих.автомобилей и многих моментов
времени привело бы к чрезвычайно громоздким
вычислениям. При решении практических задач можно допустить
некоторые упрощения ради получения .приемлемых в
вычислительном отношении алгоритмов, если при этом
ошибки -не будут превышать допустимых пределов. Мы
разделим задачу исследования автомобильного движения
на две части: 1) изучение движения одиночного
автомобиля в динамической схеме, т. е. с учетом действующих
на автомобиль сил и 2) изучение автомобильного
движения с учетом дорожного траффика, т. е. совместного
движения большого количества автомобилей в
.кинематической схеме. Другими словами, при изучении дорожного
траффика будем считать скорость автомобиля известной
функцией параметров движения. При таком подходе к
делу модель автомобильного движения существенно
упростится, так как не нужно будет на каждом шаге мо-
делирования прибегать к решению дифференциальных
уравнений, представляющих собой модель движения
одиночного автомобиля. Однако построение зависимости
скорости автомобиля от параметров движения,
необходимой для реализации кинематической схемы, можно
сделать только на основе моделирования движения
одиночного автомобиля в динамической схеме. Существуют,
правда, и экспериментальные методы построения
требуемых зависимостей.
Мы будем предполагать, что скорость свободного
движения автомобиля выражается известной функцией
1> = 1ф\/>Д/,ц,Ст) E.10)
для каждого типа автомобилей. Здесь Ст—лолный вес
автомобиля с грузом.
Соотношения для определения .скорости из формулы
E.10) встречаются »в специальной литературе.
Движение автомобилей на двухполосной дороге
представляется следующим образом. В лункте А (а также в
пункте В) в случайные моменты времени появляются
автомобили. Если дорога свободна, то первый
автомобиль пройдет весь участок АВ со скоростью свободного
движения, даже если его обгонят другие автомобили.
Если скорость второго автомобиля не больше скорости
первого, то он также пройдет весь путь со скоростью сво-
93

бодаого движения и т. д. Если скорости данного
автомобиля больше, чем скорость впереди идущего автомобиля,
и данный автомобиль успеет догна/ь .впереди идущий
(время движения до конца маршрута данного
автомобиля меньше, чем у впереди идущего), появляется
необходимость обгона. О/бгон считается возможным, если:
1) отсутствует знак, запрещающий обгон; 2) имеется
место для обгоняющего автомобиля ib своем ряду
движения длиной не менее 2/о; 3) во встречном ряду имеется
интервал не менее Lo, обеспечивающий безопасность
обгона; 4) ,во встречном ряду нет таких обгоняющих
автомобилей, которые могут помешать обгону данного
автомобиля.
Интервал Lo зависит от скоростей обгоняющего
автомобиля и скоростей движения в прямом и встречном
потоках. Если обгон невозможен, автомобиль снижает
скорость до скорости впереди идущего автомобиля и
продолжает движение ic этой скоростью.
Заметим, что обгон может произойти лишь в том
случае, если место для обгоняющего автомобиля в своем
ряду отдалено от обгоняющего автомобиля не более чем
на т автомобилей, где т определяется условиями
видимости.
Если обгон возможен, обгоняющий автомобиль
выезжает из 'Своего ряда и, двигаясь параллельно обгоняемым
автомобилям, занимает выбранное место. Далее
движение происходит аналогичным образом, пока автомобиль
не достигнет «конечной точки маршрута. Аналогично же
выглядит картина движения во встречном потоке.
20. Структура моделирующего алгоритма
Мы уже отмечали, что процесс автомобильного
движения существенно отличается от традиционных
процессов массового обслуживания. Поэтому и структура
моделирующего алгоритма обладает особенностями по
сравнению с алгоритмами, рассмотренными в предыдущей
главе.
Описываемый ниже моделирующий алгоритм состоит
из следующих основных блоков (подалгоритмов):
формирование потока автомобилей;
94

определение скорости автомобиля и его перемещения;
отыскание мест\ для обгоняющего автомобиля в своем
ряду;
оценка интервала ео 'встречном потоке;
переход к .встречному потоку;
фиксация состояний системы (результатов
моделирования).
Для представления моделирующего алгоритма в виде
операторной схемы нам потребуются следующие
операторы:
Pi — проверка наличия в регистре данных об
автомобилях (проверка усло(вия v>0, где v — количество
автомобилей, помещенных в регистр);
Ф2 — формирование параметров нового автомобиля,
появившегося на дороге (формирование момента
появления автомобиля с помощью случайных чисел, выбор
по жребию типа автомобиля, его скорости и других
характеристик);
Рз—проверка принадлежности tj — момента
появления автомобиля на дороге, данному Д^;
А4 — запись в регистр данных о сформированном
автомобиле; подсчет числа автомобилей в регистре, v+1;
А5 — выборка из регистра данных об автомобиле;
Кб —счетчик числа автомобилей в регистре
(реализует операцию v—1);
д7 — отыскание свободных ячеек для запоминания
расстояния sA, скорости автомобиля vmax, конкретной
скорости движения vv и других параметров автомобиля;
Kg — счетчик количества автомобилей, движущихся
по дороге в данном направлении (реализует операцию
N+\)\
F9 — формирование п = 0 (п — количество
автомобилей, передвижение которых в данном интервале А/ уже
рассмотрено);
Кю — счетчик количества рассмотренных
автомобилей (реализует операцию /г+1);
Рп—проверка условия n^N;
F12 — выбор номера очередного автомобиля и его
характеристик (просмотр автомобилей начинается от
пункта В для прямого потока);
А13 — определение характеристик дороги в точке
местонахождения рассматриваемого автомобиля (уклоны,
повороты, знаки -и т. д.);
95

Аи — определение vv — реальной скорости движения
автомобиля с учетом уклонов, лаворо/ов, дорожных
знаков и т. д.;
Ai5 — вычисление перемещения ^s = ap«A/ очередного
автомобиля;
Aie — вычисление расстояния sA = s^p +As для
рассматриваемого автомобиля;
Pi? — проверка условия $А<АВ;
Pis— .вычисление ,и запоминание «момента Тв
прибытия автомобиля в пункт В (для прямого потока);
Pig — проверка условия ш>0 (признак со = 1
характеризует обгоняющий автомобиль, со = 0— необгоняющий);
Р20 — проверка условия sA>sm где sM — координата
места, которое автомобиль займет после окончания
обгона;
F21 — формирование ю=0;
Р22 — проверка условия sAn^SAn-\—/б",
А2з — запись sAn в счетчик sA;
. SAn—SAn-\ + ei
А24 — вычисление б/= >'
Up
А25 — определение момента 'времени U—Ы\
А2б — определение величины сдвига sA—vv*6f = sA\
А27—сдвиг автомобилей на vV'8t;
р28 — проверка наличия знака, запрещающего обгон;
Р2э — проверка условия /е>0 (все ли впереди идущие
автомобили уже просмотрены); смысл величины k
выясняется в определении оператора Р54;
Рзо — проверка наличия расстояния, равного 2U,
между впереди идущими автомобилями;
Asi—запоминание номера автомобиля, впереди
которого разместится обгоняющий;
А32 — переход к следующему автомобилю (реализует
операцию k—1);
F33 — сн-ижение скорости автомобиля до скорости
впереди идущего;
А34 — определение sAn с учетом снижения скорости ир;
А35 — .возврат автомобилей и соответствии с ранее
вычисленными sA;
Рзс — проверка условия Т[п)<Т^~~1\
А37 — определение 6^* — момента, когда я-й и n+1-й
автомобили будут находиться на расстоянии h Друг от
друга;
96

F38 — снижение-скорости автомобиля до скорости
впереди идущего (аналогичен оператору F3s);
AG9 — вычисление перемещения автомобиля за б/*
(аналогичен оператору Л is);
А4о — вычисление расстояния 5л (аи а логичен
оператору Л16);
A4i — определение момента 7в прибытия автомобиля
в пункт В;
А42 — 'приближенная оценка номера автомобиля
встречного потока для окрестности места обгона;
А4з — обращение к встречному .потоку;
А44 — определение sA автомобиля встречного потока к
моменту t—6^;
Р45 — проверка неравенства s^CTp ^s^6r +L0
(превышает ли интервал между обгоняющим и ближайшим
встречным автомобилем величину Lo);
F46 — формирует признак f$ = l; (|3 = 0, если цикл
отбора автомобиля встречного потока состоит из одной
проверки; р=>1 в противном случае);
А47 — переход к (д+1)-му автомобилю встречного
потока;
Р48 — проверка неравенства |3>0;
А49 — переход к (п—1)-му автомобилю встречного
потока;
А5о — .вычисление интервала AsjCTp между
предыдущим и последующим автомобилями во встречном потоке;
<P5i — проверка условия As^CTp <L0;
А52—переход к (/г+1)-му автомобилю встречного
потока (аналогичен оператору Л47);
Кбз — счетчик количества п просмотренных
автомобилей (реализует операцию /г+1);
Р54 — проверка условия n^k, где k — количество
автомобилей встречного потока, находящихся в поле
видимости водителя;
Рб5 — проверка условия со>О;
Абб — определение координаты sA и величины Lo для
обгоняющего автомобиля встречного потока;
Р57 — проверка условия sA^s^v +LOj состоящая в
том, что интервал между обгоняющим автомобилем и
ближайшим автомобилем встречного потока превышает Lo;
F58 — 'сообщение автомобилю его максимальной
скорости (скорости движения без помех);
F59 — сообщение автомобилю признака со = 1;
7. Заказ 3012. 97

Або — определение sa обгоняющего автомобиля для
конечного момента интервала At;
Aei — переход к встречному потоку;
Аб2 — реализует операцию р+'1; (р формальный
параметр, используемый для перехода к встречному
потоку; при / = 0 начальное значение р = 0);
Рбз — (проверка условия р>2;
Ff,4 — формирование р = 1;
Рбз — /проверка условия р>1;
Абб — ввод начальных условий для моделирования
очередной реализации;
А67 — переход к следующему At (реализует операцию
U = /<_,+ At);
Рб8 — (проверка условия /*<Г, где Т — 'верхняя
граница интервала времени, в течение которого моделируется
процесс;
А69 — фиксация результатов моделирования за
реализацию;
Р7о — проверка условия М<М*, где М — текущее
число реализаций процесса, а М* — заданное;
Ап—обработка результатов моделирования;
Я72 — выдача результатов.
Располагая введенными операторами, запишем
операторную схему моделирующего алгоритма:
68 лК 1,8 ^к тл а9 1А jr 3,6 Л Тг2 4г? 9,23,36,41,г
Pf Ф2 Р.у|7 А4 А5 Кб А7 К8 F9 Кю
ОС 1 
Plim Fi2 А13 А14 А15 Ai6 Pjy9 Ai8 P.19122 P.20123
,-, 19,21 o , 20,22,35 Л10 22 * Л л л л t33
Г21 1 22Y24 А23 /\24 ^25 А26 А27 A^g
"Ag F» A34
ц Азт FM А39 А4о А41? 31А42 А43 43'47>49А44
F46A447445Plf ^"AsoPtf 51>55А52 К53 ^
Р 55J52 А56 Рб7|33 '' F53 F59 АбО Аб1 Кб2 Р 63J65 F64
63,64nt67 70 д 65,66 д 65,67 Dtl д р1б6 д
Рб5|67 А66 Аб7 Рб8 Аб9 ^70 А71 П72
Блок-схема моделирующего алгоритма представлен^
на рис. 14,

Описание функционирования алгоритма начнем с
блока формирования потока автомобилей. Работа его
протекает следующим образом. Оператор Pi проверяет
условие v>0, выполнение которого свидетельствует о
наличии в регистре автомобилей сформированных в
предыдущем интервале At. Если автомобили в регистре
имеются, управление по стрелке с индексом 1 передается
оператору А5 для выборки параметров автомобиля из
регистра; если нет—,тто стрелке с индексом 0 оператору
Ф2, который формирует параметры нового автомобиля.
Работа оператора <t>2 осуществляется так: из
совокупности случайных чисел, равномерно распределенных гв
интервале @, 1), выбирается очередное случайное число ?
и преобразуется к заданному закону распределения.
Полученное таким образом случайное число Z\ есть интервал
времени между моментами появления (/—1) -го и /-го
автомобилей. Поэтому момент появления /-го автомобиля на
дороге tj = tj-i + Zj. Затем с помощью процедуры «.выбор
по жребию» определяется тип автомобиля и из блока
памяти ЭВМ выбираются все интересующие нас
характеристики (максимальная скорость vmax, вес, габариты
и т. д.). Управление передается оператору Рз, который
(проверкой неравенства tj^ti) выясняет
принадлежность сформированного автомобиля рассматриваемому
интервалу At. Если неравенство tj^ti не выполнено
(момент появления автомобиля находится в данном At),
управление по стрелке с индексом 0 передается
оператору А7 для организации счетчика sA (расстояние
от пункта Л), и выбора ячеек памяти, куда
заносятся остальные сформированные данные. Оператор Ks
увеличивает на 1 величину TV — количество автомобилей на
дороге (движущихся в данном направлении). Если
неравенство, проверяемое оператором Р3, выполнено,
сформированный автомобиль не принадлежит
рассматриваемому At, данные о нем заносятся в регистр оператором
А4, который, кроме того, реализует операцию v+1 (v —
количество машин в регистре). На этом блок
формирования потока автомобилей свою работу заканчивает.
Перейдем <к описанию следующего бло<ка —
определение скорости автомобиля и его перемещения. Этот блок
работает циклически, просматривая по порядку все
автомобили на дороге. Б начале цикла п = 0 (оператор F9);
п—'количество автомобилей, которые в данном цикле
7* 99

61
Переход н встречному потону
62
Jl
о63'
69
Формирование р -- /
65^
>/}
66
67
Ввод начальных данных длл
новой реилизации
X
Лереход к следующему л t
X
63
Фиксация результатов
за реализации)
71
72
?
Оброботко
результатов моделирования
Выдача результатов
Формирование
параметров нового автомобиля
Запись даннь/х в ре-
гистр . V * /
Вы бор на данных из регисгра
—*¦! Организация счетчика >
Формирование п=0
/л . f
1?
/3
Выбор номера
очередного автомобилл
Определение
характеристик дороги
Определение реальной
скорости автомобиля
16
I Вычисление &$ = Vp &t
Вы числение $А=s/p+ д S I
18
Определение момента
прибытия в точку В
23
N6
Запись SAn в счетчик SA Г
Вычисление приращения St
25
26
Определение момента tL -
Определение сдвига
27
Сдвиг автомобилей но Vp
Рис. 14.

Продерна наличия зна -
/A f<o, запрещающего обгон
29
О
30
к>0\0__
Г/
31
Запоминание номера
адтомобиля
32
Переход к следующему
адтомобил/о, /<-/
33
Снижение скорости
адтомобиля
Грубая оцени а номера ад/по -
мобиля дстречноео по/пона
±
Переход н бстречному л о/пону
Определение SA к моменту t-frt
>/ дА
Формиродание /=/
±.
Переход н (п*!)-му адгомойилю),—
31
38
39
J
Определение
момента fit *
Снижение снорос-
ти абтомобиля
Вычисление
перемещения A S = Vp fit *
90
Вычисление
U
Определение момента
прибытия б точну В
Переход н (л-!)-му абгомобилю
50
Вы числение л s?crp= SAn -SAn^f
5/
Переход н (лЧ)-му адто-
мобилю дс/лречноео погона
59
Сообщение адтомобилю
максимальной снорости
60
Сообщение адтомо-
билю признана oj=/
Определение SA для
момента a t

уже рассмотрены. При переходе к очередному
автомобилю оператор Кш прибавляет '1 к числу п. Оператор Рп
проверкой условия n^N определяет, все ли автомобили
на дороге уже просмотрены; если условие, проверяемое
оператором Рп, .выполнено, управление по стрелке с
индексом 1 (Передается оператору ?Х2 для выбора номера
очередного автомобиля (просмотр автомобилей для
прямого потока начинается с автомобиля, ближайшего к
пункту В).
Оператор Ai3 извлекает из памяти машины
характеристики интересующего нас участка дороги (к ним
относятся уклоны, подъемы, повороты, знаки,
ограничивающие скорость, характеристики качества покрытия дороги
и все остальные величины, которые потребуются нам для
реализации оператора Аи — определение скорости
автомобиля). Определив с помощью оператора Аи скорость
автомобиля, вычисляем его перемещение As за время At
(оператор А15) и новую координату sa=sa+As
(оператор Ale). Оператор PJ7, проверяя условие sA<AB,
определяет, принадлежит ли местоположение автомобиля
рассматриваемому отрезку дороги АВ. Если sA<AB —
автомобиль не достиг еще пункта В (для прямого
потока), управление по стрелке с индексом 1 передается
оператору Pig. В противном случае определяется момент Гв
прибытия автомобиля в пункт В (оператор Ais). От
оператора Aie управление передается оператору Рзб, который
осуществляет проверку неравенства Т%<Твп~1. Здесь
Т% — момент прибытия п-то автомобиля в пункт В, а
Tq—момент прибытия (п—1)-го. Мы приняли
предположение, что на последнем (перед пунктом В) интервале
At обгон исключается. Такое предположение на практике
часто оказывается справедливым. Кроме того, мы не
располагаем сведениями о движении автомобилей за точкой
В. и поэтому не можем проверить возможность обгона на
этом интервале At, Вследствие сказанного
справедливость неравенства Т%<Т%~1 означает лишь формальное
изменение порядка следования (п—>1)чго и я-го
автомобилей; поэтому в некоторый момент времени U + 6t* п-й
автомобиль должен снизить скорость до скорости
впереди идущего. От Рзб по стрелке с индексом 1 управление
передается оператору А37- Оператор А37 определяет
момент 6/*, коцда данный и впереди идущий автомобили
102

будут находиться на расстоянии k (расстояние
безопасности движения). Этот момент интересен с той точки
зрения, что до 6^* п-и автомобиль будет двигаться со авоей
максимальной скоростью, а после 6/*—со скоростью
впереди идущего автомобиля. Оператор F38 присваивает
данному автомобилю скорость впереди идущего.
Оператор Азэ вычисляет новое As = fp6i* (работает аналогично
оператору Ai5). Оператор А4о вычисляет новую
координату sAn = sA + As (аналогично оператору Ai6) и, наконец,
оператор A4i вычисляет момент прибытия
рассматриваемого автомобиля в .пункт Я. От оператора А41 управление
передается оператору Кш Для подсчета количества
рассмотренных автомобилей; туда же по стрелке с
индексом 1 передается управление и от оператора *Рзб, если
неравенство Т%<Т^~1 не выполнено (я-й автомобиль
прибывает в В, не нарушая порядка движения).
Теперь вернемся к оператору Pi7. Если автомобиль
еще не достиг пункта В (sA<AB), ino стрелке с
индексом 1 управление передается оператору Pig, который
проверяет неравенство со>0. Признак со=1 присваивается
обгоняющему автомобилю, поэтому выполнение
неравенства со>О свидетельствует о том, что рассматриваемый
автомобиль является обгоняющим. Тогда по стрелке с
индексом 1 управление передается оператору Р2о,
который, проверяя нера>венство sA>sM/ определяет, достиг ли
данный обгоняющий автомобиль места, на которое он
должен стать, закончив обгон. Если неравенство sA>sM
выполнено (обгон закончен), оператор F2i присваивает
данному автомобилю признак со=0 (что соответствует
необгоняющему автомобилю) и передает управление
оператору Р22. Если обгон еще не закончен sA<sAT,
управление по стреше с индексом 0 передается оператору А2з,
который записывает в счетчик sA данного автомобиля
величину sAn, вычисленную оператором Ai6, и передает
управление оператору Кю для перехода к следующему
автомобилю. Оператор Р22 проверяет неравенство
sAn^sAn-i—/5; смысл этой проверки заключается в том,
чтобы определить, имеется ли между рассматриваемым и
впереди идущим автомобилем расстояние, равное или
большее U (расстояние безопасности). Если sAn^
^sAn-\ — lQ (т. е. рассматриваемый автомобиль в своем
движении не догнал еще предыдущего), управление по
стрелке с индексом 1 передается оператору А23 для запи-
103

си sAn в счетчик sA и затем оператору Кш для перехода
<к новому автомобилю. Если sAn>sAn-\ — h> т. е. возникла
ситуация, соответствующая ложности неравенства
Т% <Твп~1 (оператор Рзб), управление по стрелке с
индексом 0 передается оператору А24. Оператор А24 определяет
момент времени 6/, когда эти автомобили будут
находиться на расстоянии h друг от друга (dt отсчитывается
от начала данного АО- Оператор А25 вычисляет этот
момент относительно /г=0 (начала отсчета). Оператор А26
определяет местоположение автомобиля к моменту
ti — At; sA*=sA—vv*6t. Итак, если неравенство sAn^.
^.sAn-i — /б не выполнено, т. е. впереди идущий
автомобиль .мешаег данному двигаться со своей максимальной
скоростью, данный автомобиль вынужден идти на обгон.
Чтобы определить воЗ(Можность обгона, требуется решить
несколько вопросов и в первую очередь вопрос о наличии
места в своем ряду, куда обгоняющий автомобиль может
стать по окончании обгона. Для отыскания места в своем
ряду (это место ищется среди k+\ впереди идущих
автомобилей) необходимо переместить k + \ впереди идущих
автомобилей на vp-6t назад, т. е. восстановить на дороге
ситуацию, возникшую к моменту ti — Ы. Этим и
занимается оператор А27.
На ЭТО1М мы закончим рассмотрение блока
определения скорости автомобиля и его перемещения и перейдем
к описанию блока проверки возможности обгона.
Оператор Р28—первый оператор блока проверки возможности
обгона — решает этот вопрос, выясняя наличие
дорожного знака; запрещающего обгон на данном участке
дороги (обгон также запрещен на спусках, поворотах и
других особо опасных участках дороги). Если обгон
запрещен, управление по стрелке с индексом 1 передается
оператору F33 для сообщения данному автомобилю
скорости впереди идущего; если разрешен — по стрелке с
индексом 0 оператору Р2д, которым начинается цикл
отыскания места в своем ряду (поиск места ведется среди k
впереди идущих автомобилей). Оператор Р2э проверяет,
все ли впереди идущие автомобили уже просмотрены.
Если k>0 — автомобили еще имеются, управление
передается (по стрелке с индексом 1) оператору Рзо для
проверки неравенства sAk + 2k^sAk-] (проверки наличия
интервала, длиной 2/б между k-ы и k~-\-u автомобилем).
Если неравенство sAk + 2k^sAk~i не выполнено, управ-
104

ление ,по стрелке с индексом 0 лередается оператору А32
для перехода к следующему автомобилю (реализуется
операция k— 1). Если интервал длиной 2Ц не обнаружен
среди всех k впереди идущих автомобилей, это значит,
что данный автомобиль не может совершить Oi6roH. Тогда
оператор F33 сообщает ему скорость впереди идущего.
Оператор А34 определяет sA автомобиля с учетом
измененной скорости vp, а оператор А35 воз-вращает k+ 1
сдвинутые автомобили на свои прежние места, т. е.
восстанавливает на дороге ситуацию, имевшую место к моменту
t%. От оператора А35 управление передается оператору А2з
для записи sAn в счетчик sA.
Возвратимся теперь к оператору Рзо- Если
неравенство SAh+Qh^SAk-u проверяемое оператором Р3о,
выполнено, управление по стрелке с индексом 1 передается
оператору А3ь который фиксирует номер автомобиля,
впереди которого имеется место для обгоняющего
автомобиля. От оператора A3i управление передается
оператору А42, которым начинается новый блок — оценка
интервала во встречном потоке.
Оператор А42 приближенно оценивает номер
автомобиля встречного потока, ближайшего к обгоняющему.
Оператор А4з осуществляет переход к встречному потоку
для некоторой окрестности места обгона. Оператор А44
определяет sA автомобиля (номер которого был получен
в операторе А42) к моменту ti — 6t. Для этого же
автомобиля оператор Р45 проверяет неравенство s?CTp ^5дбг +
+ L0, т. е. определяет, превышает ли интервал между
обгоняющим и ближайшим встречным автомобилем
величину Lq. Если условие, проверяемое оператором Р45,
выполнено, управление по стрелке с индексом 1
передается оператору Р48. Однако остается открытым вопрос о
наличии автомобилей внутри интервала Lo. Для проверки
последнего обстоятельства придется последовательно
переходить к (п—1)-му или (п+1) встречным
автомобилям. От оператора Р45 по стрелке с индексом 0 управление
передается оператору F46 для формирования признака
р= 1, Оператор А47 обеспечивает переход к (я+ 1)-му
автомобилю встречного потока и опять передает управление
оператору А44. На этом заканчивается первый полуцикл
обследования автомобилей встречного потока. От
оператора Р48 (если условие |3>0 не выполнено) по стрелке с
индексом 0 управление передается оператору А4э для пе-
105

рехода к (п—1)-му автомобилю встречного потока, затем
управление опять передается оператору А44.
На этом заканчивается второй из упомянутых
полуциклов.
В случае, когда условие, проверяемое оператором Р48,
выполнено, управление по стрелке с индексом 1
передается оператору Aso Для вычисления фактического
интервала между автомобилями встречного штока AS^CTp =
= 5лп—5лт2-1. Оператор Psi проверяет условие,
состоящее в том, что вычисленный интервал меньше, чем Lo.
Если это условие выполнено и, следовательно, обгон
невозможен, упрлзление передается оператору Рзз для
присвоения дань ому автомобилю скорости впереди
идущего. Если А5а>^о — обгон возможен, управление по
стрелке с индексом 0 передается оператору А52, которым
начинается подалгоритм, определяющий, не мешает ли
обгону обгоняющая встречная машина. Этот
подалгоритм работает циклически, просматривая, начиная с п-й
автомашины, k автомашин встречного потока.
Оператор А52 осуществляет переход к n+il-й автомашине;
оператор К53 подсчитывает л — количество уже
просмотренных автомобилей. Оператор Р54 проверяет, все ли k
автомобилей мы уже просмотрели и, если n<k, по стрелке с
индексом 0 передает управление оператору Р55.
Оператор Рб5> проверяя значение признака о, определяет,
является ли рассматриваемая автомашина обгоняющей или
нет. Если 0 = 0 управление передается снова оператору
Arj2 для перехода к следующей автомашине. Если со>0 —
автомобиль обгоняющий, оператор Абб вычисляет 5а и Lo
и затем оператор Р57 проверкой неравенства s? ^s?6r +
+ L? , определяет, помешает ли он обгону или нет. Когда
нераьенство s? >sj?6r + L0. оказывается
невыполненным, т. е. встречный обгоняющий автомобиль не позволит
осуществить обгон, управление по стрелке с индексом
0 передается снова оператору F33, который присваивает
данному автомобилю скорость впереди идущего. Если
sa =^а6г -Ь^-о — обгон возможен, по стрелке с
индексом 1 управление передается оператору F58. Оператор F58
сообщает данному автомобилю его максимальную
скорость, т. е. ту скорость, с которой автомобиль движется,
если нет помех дзижению. Затем оператор F59
присваивает данному автомобилю признак со = 1 и оператор Або
106

определяет sAn— координату обгоняющего автомобиля
для конечного момента интервала А/. От оператора Ago
управление передается оператору А35, который,
возвращая k+l впереди идущих автомобилей вперед,
восстанавливает нормальную ситуацию на дороге. Затем
управление передается оператору А2з для записи
вычисленного оператором А605ап в счетчик sA и затем оператору
Кю Для перехода «к новому автомобилю.
На этом заканчивается описание работы алгоритма
над машинами прямого потока (из пункта А в пункт В).
Оставшаяся группа операторов занимается переходом к
встречному потоку, переходом к новому At, фиксацией и
обработкой результатов моделирования. Рассмотрим
подообнее работу этого последнего блока.
Оператор Аб1 обеопеч'ивает переход к встречному
потоку и .передает управление оператору Кб2, который
прибавляет единицу к параметру р. Оператор Рбз проверяет
условие р>2. Когда это условие выполнено, управление
по стрелке с индексом 1 передается оператору F64,
который формирует р=1 и передает управление оператору
Рб5- Если условие, проверяемое оператором Р6з, не
выполнено, по стрелке с индексом 0 управление
непосредственно передается оператору Рб5- Этот оператор
проверяет условие р>1; если условие выполнено, по стрелке с
индексом 1 управление передается оператору А67. В
случае, если условие, проверяемое оператором Рб5, не
выполнено, по стрелке с индексом 0 управление передается
оператору Р68.
Смысл этих процедур состоит в следующем.
Предположим сначала, что оператор Aei вступает в действие
впервые за реализацию. В этом случае р=0 (начальное
условие); тогда от оператора Р6з мы попадаем к
операторам Рб5 и Р68, т. е. просматриваем тот же интервал А^
для встречного потока. Если же речь идет о возвращении
к прямому потоку, то р>1 (оператор Кб?) и от
операторов Рбз и Р65 мы попадаем к оператору А67, т. е.
переходим к очередному At. При следующем переходе к
встречному потоку р>2, но это исправляется оператором
F64. Оператор А67 обеспечивает переход к следующему
At и передает управление оператору Рее, который лро;ве-
ряет условие ti<T. Если это условие не выполнено,
значит интервал моделирования [О, Т] закончен и тогда по
стрелке с индексом 0 мы переходим к оператору А6о,
107

фиксирующему результаты моделирования за
реализацию. Если же условие, проверяемое оператором Рб8,
выполнено, мы переходим к моделированию процесса и уп-
равлеиие передается оператору Рь
От оператора А6о управление передается оператору
Р7о. Если условие, проверяемое оператором Р7о,
выполнено, управление передается по стрелке с индексом 1
оператору Абб для ввода начальных условий,
соответствующих переходу к новой реализации процесса. Далее
управление передается упомянутому ранее оператору АгO.
В случае, если условие, проверяемое оператором Р7о, не
выполнено, работает стандартная пара операторов A7i и
Я72 — обработка и выдача результатов.
В заключение остановимся кратко па способах
приближенной оценки номера автомобиля встречного потока
(оператор А42). Для проверки возможности обгона
необходимо оценить расстояния между соседними
автомобилями встречного потока, ближайшими к автомобилю,
претендующему на обгон. Если для этой цели
пользоваться перебором Есех автомобилей, мы придем к
непомерно большому количеству операций на ЭВМ.
Желательно значительно сократить количество
рассматриваемых при переборе автомобилей. Для этого необходимо
хотя бы примерно оценить номер автомобиля встречного
потока, находящегося поблизости от автомобиля,
претендующего на обгон. Здесь могут быть пригодны
различные способы, в зависимости от особенностей модели,
В частности, в рассматриваемом моделирующем
алгоритме применялся следующий прием.
Пусть s' — координата автомобиля, претендующего
на обгон в момент времени t\ Расстояние АВ—s'
встречный автомобиль, движущийся со скоростью уСр, пройдет
за интервал временя
rAB-s'
т= •
Если ?,встр — интенсивность встречного потока
(автомобилен в час), то за интервал 'времени f—?0—т, где /о—
момент начала движения автомобилей на дороге, через
точку s' пройдет приблизительно
П = Х13стР(//-/о~т)
автомобилей. Целая часть числа п является
приближенной оценкой номера автомобиля встречного потока,
ближайшего к автомобилю, претендующему на обгон.
108

Заметим, что разработка более точных способов
оценки числа п способствовала бы сокращению
непроизводительного расхода операций ЭВМ при моделировании.
21. Исходные данные и результаты моделирования
Первичные исходные данные для моделирования
процесса автомобильного движения могут быть
непосредственно перечислены исходя из рассмотрения
операторов моделирующего алгоритма и всех входящих в них
расчетных формул. К ним относятся:
характеристики дороги (длины участков между
особыми точками, величины подъемов и уклонов,
координаты и содержание дорожных знаков, значения
коэффициентов сцепления и поправок к ним и т. д.);
длительности интервала моделирования (/о, Т) и
интервалов А/;
законы распределения моментов -появления
автомобилей на дороге;
"количестзо типов автомобилей и состав потока по
типам (процент);
законы распределения параметров автомобилей (вес,
вес груза, максимальная скорость и т. д.) по типам
автомобилей;
параметры, входящие в выражения для ^величии Iq, Lo,
vp, к и т. д.
В результате моделирования на печать обычно
выдаются:
средняя длительность движения автомобилей данного
типа от пункта А до лункта В;
средняя потеря времени движения автомобилей
данного типа на маршруте АВ по сравнению с временем
свободного движения;
средняя реальная скорость ир автомобилей данного
типа на маршруте АВ;
относительный коэффициент av влияния
интенсивности
Vp
для автомобилей данного типа; здесь прасч — средняя
реальная скорость при интенсивности движения, близкой
к нулю;
109

среднее количество несостоявшихся обгонов за
реализацию процесса по типам автомобилей (иногда по
причинам: из-за дорожных знаков, отсутствия места в своем
ряду, интервала во встречном потоке и т. д.);
среднее количество состоявшихся обгонов по типам
автомобилей и ряд других.
Рассмотренная модель может быть пригодной для
решения многих практических задач, возникающих при
проектировании автомобилей и дорог, а также при
организации перевозок. Наиболее распространенной является
задача построения зависимости (например, графика)
коэффициента av от интенсивности движения X. Эта
задача может быть решена путем моделирования процесса
при заданных значениях X и фиксированных значениях
других (параметров (состав потока, характеристики
автомобилей и т. д.). Обычно для практики представляет
интерес набор упомянутых (Графиков, соответствующий
различным значениям состава потока и характеристик
автомобилей. Иногда также строятся наборы графиков для
типовых дорот с различными параметрами.
Поскольку наборы таких графиков, охватывающие
значения параметров в пределах, достаточных для
решения практических задач, были бы слишком громоздкими
(или, в противном случае давали бы лишь
приближенную оценку зависимости av(X) для различных условий
движения), целесообразно использовать модель для
получения зависимости av(X) 'при конкретных условиях,
представляющих интерес для какой-нибудь .практической
задачи. Тогда модель была бы инструментом
использования в организациях, занимающихся проектированием
автомобилей, дорог или организацией перевозок.
Помимо оценки зависимости av(X) три
фиксированных параметрах, существенный интерес представляет
исследование ее при варьированных параметрах. Обычно
таким путем удается выяснить степень стабильности
зависимости av(X)9 а также изменение ее характера при
изменении отдельных параметров.
Сведения о состояниях (процесса, получаемые методом
моделирования, .позволяют решать и другие задачи,
представляющие (Практический интерес. Здесь можно отметить
задачу об .изучении переходного режима, связанного с
постепенным насыщением дороги автомобилями.
Решение этой задачи позволяет оценить аффективные значе-

ния интенсивности движения, соответствующие заданным
условиям, которые более точно согласуются с опытом,
чем формальная зависимость av(X).
Весьма интересные (хотя, быть может, более
трудные) задачи возникают при так называемом системном
подходе к изучению процесса автомобильного движения.
Рассматривая этот «процесс как систему, можно ставить
вопросы об оценке его эффективности (с различных точек
зрения), надежности (с учетом отказов автомобилей и
выходов из строя элементов дороги),
помехозащищенности (с учетом действия погоды, колонн медленно
движущихся машин и других помех), а также устойчивости по
отношению к различным возмущениям (увеличение
интенсивности, снижение скорости автомобилей и
проходимости отдельных участков дороги и т. д.).
Литература
¦1. Проблемы научной организации управления
социалистической промышленностью. Сб. статей под ред. Д. М. Гвишиани и др.
М., «Экономика», 1968.
2. Г. А. А л и е в* и др. Моделирование производственного
процесса автоматизированного стана печной сварки труб. «Проблемы
кибернетики». Вып. 9, М., Физматгиз, 1963.
3. Н. П. Б у с л е н к о. Математическое моделирование
производственных процессов на цифровых вычислительных машинах. М.,
«Наука», 1964.
4. Н. П. Б у с л е н к о. Моделирование сложных систем. М.,
«Наука», 1968.
5. Е. С. В е и т ц е л ь. Теория вероятностей. Издание 4. М.,
«Наука», -1969.
6. Е. С. В е н т ц е л ь. Введение в исследование операций. М.,
«Советское радио», 1964.
7. И. Н. Коваленко. О некоторых классах сложных систем.
Известия АН СССР. «Техническая кибернетика», № б, 1964; № 1,
3, 1965.
8. И. А. П л я ц и д е в с к и й. Информационные системы в
технике и экономике. М., «Московский рабочий», 1966.
9. И. В. Д у н и н-Б а р ко в ски й, Н. В. Смирнов. Теория
вероятностей и математическая статистика в технике (общая часть).
Гостехиздат, 1955.
10. И. Г. Венеции й, Г. С. Кильдишев. Основы теории
вероятностей и математической статистики. М., «Статистика», 1968.
11. А. Я. Боярский. Математика для экономистов. Издание
2. М., «Статистика», 1961.

ОГЛАВЛЕНИЕ
Введение 3
Глава I. Метод статистических испытаний (Моите Карло) 7
1. Сущность метода 7
2. Попадание случайной величины в заданный интервал.
Среднее значение функции от случайной величины . 14
3. Кратные -интегралы 20
4. Требуемое количество операций 22
Глава II. Получение и цреобразовашие случайных чисел . . 27
5. Случайные числа 27
6. Получение случайных чисел с заданным законом
•распределения 31
7. Приближенные способы преобразования случайных
чисел , 36
8. Моделирование случайных векторов 39
Глава III. Применение метода статистических испытаний
для моделирования сложных систем 44
'9. Математическая модель 44
10. Этапы формализации реальных сложных систем . 47
11. Операторная схема моделирующего алгоритма . . 50
12. Метод статистического моделирования как аппарат
анализа сложных систем ....... л ... 54
Глава IV. Статистические модели систем массового
обслуживания t 58
13. Понятие системы массового обслуживания . , . „ 58
14. Моделирование потоков заявок 64
15. Моделирование одноканальной системы массового
обслуживания 68
10. Моделирование многоканальной -системы 73
17. Модификация систем массового обслуживания ... 77
Глава V. Моделирование автомобильного движения .... 87
18. Вводные замечания. Формальное описание дороги . 87
19. Описание потока автомобилей 90
20. Структура моделирующего алгоритма 94
21. Исходные данные и результаты моделирования . . . 109
Литература ь , 111
Николай Пантелеймопович Бусленко
Метод статистического моделирования
Научный редактор Н. И. Щедрин
Редактор В. Н. Третьякова
Техн. редактор Т. Д. Алексеева. Корректор Н. С. Хорошилвва
Обложка художника Л. С. Эрмаиа
Худ. редактор Т. В. Стихно
Сдано в набор 29/IV 1970 г. Подписано к печати 24/IX 1970 г. Формат бумаги
84х108/й- Бумага № 2. Объем 3,5 печ. л. Уч-изд. л. 5,63. Тираж 9700 экз.
АН 112. (Тематич. план 1970 г. № 22).
Издательство «Статистика», Москва, ул. Кирова, 39
Заказ №-3012. Областная типография Ивановского управления по печати,
г. Иваново, Типографская, 6. 11ена 34 коп.