Статистическое моделирование

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

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

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

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

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

3.1. Сущность имитационного моделирования

Сущность имитационного моделирования рассмотрим на примере.

Пример 3.1. По объекту наносится одиночный ракетный удар. Радиус поражения $$R$$.

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

по дальности $$\sigma _{x}$$ ;

по направлению $$\sigma_y $$.

Цель будет уничтожена, если расстояние $$r $$ от нее (то есть от точки прицеливания) до центра взрыва ракеты будет меньше или равно $$R$$, то есть $$r\le R$$.Так как $$R \gg $$ размеров объекта, то цель можно считать точечной.

Наличие рассеивания исключает однозначный ответ: "цель поражена - цель не поражена". Задача носит вероятностный характер, поэтому в результате моделирования может быть получен ответ: цель будет поражена с вероятностью $$P $$.

Цель моделирования: определить вероятность $$P = P(r \le R )$$ поражения объекта одиночным ракетным ударом.

Решение

Построим декартову систему координат так, чтобы точечный объект находился в начале координат, а направление пуска ракеты совпадало с осью $$x $$ (рис. 3.1).

Возьмем две последовательности нормально распределенных случайных чисел:

$$x _{1}, x_{2} , \ldots , x_{i}, \ldots , x_{N};\\ y_{1},y_{2},\ldots ,y_{i},\ldots ,y_{N}$$ (рис 3.1) Иллюстрация к нанесению удара

Первая последовательность соответствует распределению $$М[х] = 0, \sigma_{Х} = a$$, вторая - $$М[у] = 0, \sigma_y = b$$. Матожидания $$M[x]$$,

$$М[у] $$ взяты равными нулю, так как объект поражения (точка прицеливания) находится в начале координат, то есть имеет координаты $$х = 0$$ и $$у = 0$$.

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

Моделирование

Имитируем удар, то есть мысленно нанесем удар по объекту путем определения координат взрыва. В силу идентичности закона рассеивания и его характеристик с законами распределения случайных чисел такими координатами могут быть $$x_{1}$$ и $$у_{1}$$, взятые из последовательностей случайных чисел.

Вычислим расстояние $$r _{1}$$ от места взрыва ракеты до цели:$$r_1 =\sqrt{x_1^2+y_1^2}$$

Оценим результаты имитации удара, то есть установим факт поражения или непоражения объекта:

если $$r_{1} \le R$$, то объект поражен;

если $$r _{1} > R $$, то объект непоражен.

Если объект поражен, запомним этот факт увеличением $$М $$ на единицу, то есть $$М =М +1$$ (в начале $$М = 0$$ ).

Для нахождения вероятности поражения объекта повторим имитацию нанесения удара $$N $$ раз.

Оценим вероятность через частость поражения объекта: $$P(r\le R) = \cfrac{M}{N}$$

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

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

В рассмотренном примере 3.1 при $$R = 2км$$, $$\sigma _{x} =1,5 км$$, $$\sigma = 0,8км$$ оценки вероятностей $$P $$ поражения цели при различном числе реализаций модели показаны в табл. 3.1.

Оценки вероятностей поражения цели

$$N$$	20	200	2000	10000
$$P$$	0,8	0,75	0,7615	0,7644

При $$N = 10000$$, тех же характеристиках рассеивания и других радиусах поражения $$R $$ получим:

$$R = 1\to P = 0,33; R = 1,5 \to P = 0,576; R = 2,5 \to P = 0,8823.$$

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

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

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

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

Вычисляются параметры исхода процесса при случайных значениях характеристик, полученных на этапе 2, и запоминаются. Этапы 2 и 3 соответствуют одному статистическому испытанию.

В результате $$N $$ статистических испытаний (повторений этапов 2 и 3) получают $$N $$ значений параметров исхода процесса. Вероятностные характеристики параметров исхода процесса получают в результате статистической обработки полученных случайных величин.

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

Рассмотрим еще один пример.

Пример 3.2. Транспорт 1 с грузом отправился из пункта А в пункт С через пункт В. Одновременно из пункта D в пункт Е через пункт В отправился транспорт 2. Скорости движения транспортов распределены по нормальному закону с математическими ожиданиями $$V_{1}$$ и $$V_{2}$$ и стандартными отклонениями $$\sigma_{1}$$ и $$\sigma_{2}.$$

Построить алгоритм имитационной модели (ИМ) с целью определения вероятности встречи транспортов 1 и 2 в пункте В. Расстояние от пункта А до пункта В $$S_{1},$$ а от пункта D до пункта В - $$S_{2}$$.Событие встречи считать состоявшимся, если их времена прибытия в пункт В либо равны, либо отличаются на величину, не превышающую $$\Delta t $$.

Решение

Построим схему движения транспортов 1 и 2 (рис. 3.2).

Возьмем две последовательности нормально распределенных случайных чисел:

$$V_{11}, V_{12} , \ldots , V_{1i} _{}, \ldots , V_{1N};\\ V_{21}, V_{22} , \ldots , V_{2i} _{}, \ldots , V_{2N},$$

характеристики которых соответствуют матожиданиям и стандартным отклонениям скоростей движения транспортов 1 и 2.

(рис 3.2) Схема движения транспортов

Имитируем движение транспортов 1 и 2 до пункта В со скоростями $$V_{11}$$ и $$V_{21}$$ соответственно, взятыми из последовательностей нормально распределенных случайных чисел.

Вычислим время $$t_{1}$$ и $$t_{2}$$ прибытия в пункт В транспортов 1 и 2 соответственно:$$t_1=\cfrac{S_1}{V_{11}},\,t_2=\cfrac{S_2}{V_{21}}$$

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

если $$|t_{1} -t_{2}| \le \Delta t$$, встреча состоялась;

если $$|t_{1} -t_{2}| > \Delta t$$, встреча не состоялась.

Если встреча состоялась, зафиксируем этот факт увеличением значения $$М $$ на $$1$$, т. е. $$M=M + 1 $$ (вначале $$M = 0$$ ).

Для нахождения вероятности встречи транспортов 1 и 2 повторим имитацию их движения $$N $$ раз.

Рассчитаем вероятность встречи:

$$P(\Delta t\le |t_1 - t_2|) =\cfrac{M}{N}$$

Результаты моделирования при $$N = 10 000$$ и характеристиках движения транспортов: $$V_{1} = V_{2} = 25 км/час$$, $$\sigma_{1} = \sigma_{2} = 3 км/час$$, $$S1 = S_{2} = 50 км/час$$:

$$\Delta t = 0,2час \to P =0,0908;\\ \Delta t = 0,3час \to P =0,1548;\\ \Delta t = 0,4час \to P =0,2543;\\ \Delta t = 0,5час \to P =0,3779;\\ \Delta t = 0,6час \to P =0,5571.$$

Очевидно, изложенный процесс имитации легко может быть реализован на компьютере. Представим алгоритмы моделей примеров 3.1 и 3.2 схемами (рис. 3.3 и 3.4).

В рассмотренных примерах исследуются различные процессы. Но алгоритмы моделей этих процессов (для сравнения рядом с алгоритмом задачи 3.2 (рис. 3.3) показан и алгоритм задачи 3.1 (рис. 3.4) имеют общую, практически идентичную часть (блоки 1, 5…8, на рис. 3.3 и 3.4 они выделены) и часть, которая непосредственно имитирует исследуемый процесс (блоки 2… 4).

(рис 3.3)

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

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

Алгоритм ИМ приведен на рис. 3.5. На нем выделены блоки 1, 8…11, выполняющие те же функции, что блоки 1, 5…8 в алгоритмах ИМ на рис. 3.3 и 3.4. Блоки 2…7 непосредственно имитируют нанесение удара по объекту, т. е. выполняют одну реализацию (один прогон модели). В блоке 2 переменной $$k$$ присваивается начальное число пусков ракет. Далее эта переменная используется для организации внутреннего цикла по числу пусков. После каждого пуска значение k уменьшается на 1 (блок 7). При $$k = 0$$ (блок 3) реализация завершается. Завершается она также и при поражении объекта (блок 6). Но при этом предварительно значение переменной $$М $$ увеличивается на $$1$$. По завершении $$N_{0}$$ реализаций рассчитывается оценка математического ожидания вероятности поражения объекта тремя последовательными пусками ракет.

(рис 3.5) Алгоритм модели нанесения удара тремя ракетами

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

Иногда статистическое моделирование называют "метод Монте-Карло", по городу, где процветает игра в рулетку, исход которой случаен и образуется своеобразным датчиком случайных исходов - рулеткой.

3.2. Общая характеристика метода имитационного моделирования

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

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

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

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

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

Если сложность аналитической модели с усложнением моделируемого объекта возрастает с ускорением, как показано на рис. 3.6, то сложность имитационной модели, начиная с некоторого уровня $$S_{0},$$ растет незначительно.

(рис 3.6) Иллюстрация роста сложности моделей

К достоинствам имитационных моделей можно также отнести:

простота алгоритма;

малая связность алгоритма;

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

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

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

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

3.3. Статистическое моделирование при решении детерминированных задач

Метод статистических испытаний может быть использован как численный метод решения математических задач. Именно в таком качестве он был применен в США в 1944 г. Джоном фон Нейманом при расчетах по созданию ядерного реактора.

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

Пример 3.4. Пусть $$0\le f(x)\le 1$$, $$0\le x\le 1$$. Полагаем, что функция $$f(x) $$ такова, что интеграл относится к "неберущимся".

Требуется вычислить $$S = \int\limits_{0}^{1} f(x)dx$$.

Решение

Представим функцию в координатах $$f(x) $$ и $$x $$ как показано на рис. 3.7. Как известно, численное значение интеграла данного вида равно площади $$S$$. Площадь $$S $$ состоит из множества элементарных площадок - точек. Количество точек в этой площади и будет численным значением искомого интеграла.

Имитируем координаты каждой точки значениями $$x_{i} $$ и $$y_{i}$$, принадлежащими равномерному распределению на участке $$[0, 1]$$:

$$x _{1}, x_{2} , \ldots , x_{i}, \ldots , x_{N};\\ y_{1},y_{2},\ldots ,y_{i},\ldots ,y_{N}_{} $$(рис 3.7) Вычисление интеграла

Рассмотрим пару чисел $$x_{i}, y_{i}$$. Вычислим $$f(x_{i})$$ и сравним с $$y_{i}$$. Если $$f(x_{i}) < y_{i}$$, то это означает, что точка $$M(x_{i}, y_{i})$$ принадлежит площади $$S$$. Если $$f(x_{i} )>y_{i}$$, то это означает, что точка $$M(x_{i}, y_{i})$$ не принадлежит площади $$S $$.

Введем:

$$z_i = \left \{ \begin{array}{l} 1, f(x_i)\le y_i,\\ 0, f(x_i)> y_i, \end{array}$$

Число точек, попавших в границы $$S $$ равно $$\sum\limits_{i=1}^{N}z_{i} $$, где $$N $$ - общее число точек, попавших в единичную площадь существования функции и аргумента. Отсюда следует:

$$\cfrac{\sum\limits_{i=1}^{N}z_i}{N}\approx S = \int\limits_0^1{f(x)dx}.$$

Чем больше будет элементарных площадей - точек, тем точнее будет вычислен интеграл. Приведенное решение примера справедливо для единичных областей существования функции и аргумента. Однако это несущественно, так как произвольные границы существования $$a \le x\le b,0\le f(x)\le A$$ заменой переменных можно свести к единичным границам.

Известны статистические алгоритмы численного решения многократных интегралов.

Пример 3.5. Найти оценку $$I*$$ интеграла $$I = \int\limits_0^1{dx\int\limits_x^1{(x + y)dy}} $$.

Решение

Область интегрирования ограничена линиями $$y = x$$, $$y = 1$$, $$x = 0 $$, т. е. принадлежит единичному квадрату (рис. 3.8).

(рис 3.8) Иллюстрация к примеру 3.5

Площадь области интегрирования (прямоугольного треугольника) $$S = (1 ? 1) / 2 = 0,5.$$ Используем формулу

$$I^*=S\cdot\cfrac{\sum\limits_{i=1}^{M}{f(x_i,y_i)}}{M}= 0.5\cdot\cfrac{\sum\limits_{i=1}^{M}{f(x_i,y_i)}}{M},$$

в которой $$M $$ - число случайных точек $$(x_{i}, y_{i}) $$, принадлежащих области интегрирования. У этих точек $$y_{i} \ge x_{i}$$. Если данное условие выполняется, то вычисляется

$$\sum\limits_{i=1}^{M}{f(x_i,y_i)}= {\sum\limits_{i=1}^{M}{f(x_i,y_i)}} + (x_i + y_i)$$

а число случайных точек $$M $$ увеличивается на $$1$$: $$M = M +1$$.

Результаты моделирования приведены в табл. 3.2.

Из данных табл. 3.2 (верхние пять строк) видно, что с увеличением числа реализаций $$N $$ ошибка $$\Delta S = S-I^* $$ в определении оценки интеграла $$I^*$$ уменьшается и при $$N = 1000 000$$ становится равной нулю.

Результаты моделирования примера 3.5

$$N$$	10	1000	10000	100000	1000000
$$M$$	5	500	5038	49658	500364
$$\ sum\limits_{i=1}^{M}{f(x_i,y_i)}$$	4,773	487,695	5006,152	49533,242	500191,650
$$I^*$$	0,477	0,488	0,497	0,499	0,500
$$\Delta S = S-I^*$$	0,023	0,012	0,003	0,001	0
$$M$$	6	503	4935	49833
$$\ sum\limits_{i=1}^{M}{f(x_i,y_i)}$$	5,025	494,593	4917,236	49802,019
$$I^*$$	0,419	0,492	0,498	0,500
$$\Delta S = S-I^*$$	0,081	0,008	0,002	0

В четырех нижних строках табл. 3.2 приведены результаты моделирования с другими начальными числами генераторов равномерно распределенных случайных чисел. Как видно, ошибка в оценке интеграла равна нулю уже при $$N= 100 000$$ реализаций модели.

В заключение отметим, что имитационное (статистическое) моделирование целесообразно применять в случаях:

когда нет законченной математической постановки задачи;

когда нет аналитических методов решения сформулированной задачи;

когда аналитические методы есть, но они не удовлетворяют требованиям точности и достоверности;

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

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

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

3.4. Моделирование равномерно распределенной случайной величины

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

Определение. Непрерывная случайная величина $$\gamma $$ имеет равномерное распределение в интервале $$[a, b]$$, если ее плотность вероятности $$f(x) $$ определяется так (рис. 3.9):

(рис 3.9) Плотность вероятности равномерного распределения

Значения характеристик равномерного закона распределения:

математическое ожидание $$M [ \gamma ] =\cfrac{a+b}{2}$$ ;

дисперсия $$D[ \gamma] =\cfrac{(b-a)^2}{12}$$.

При моделировании часто используются случайные числа из интервала $$[0,1]$$. Непрерывная случайная величина $$\gamma $$ равномерно распределена в интервале $$[0,1]$$, если:

$$f(x) =\left \{\begin{array}{l} 1,0\le x\le 1,\\ 0,x\prec 0,x\succ 1 \end{array}$$

В этом случае $$M[ \gamma ] =\cfrac{1}{2},D[ \gamma ] =\cfrac{1}{12}$$.

Случайное число $$x_{i} $$ из интервала $$[0,1]$$ легко преобразуется в случайное число $$x_{i}^{'}$$ для интервала $$[a, b]$$:

$$x_{i}^{?} = (b - a)\cdot x_{i} + a .$$

Применительно к двоичным дробям случайное число из интервала $$[0,1]$$ представляет собой бесконечную дробь:

$$\gamma=z_1\cdot 2^{-1}+z_2\cdot 2^{-2}+\ldots z_i\cdot 2^{-i}+\ldots;\,\,i=\overline{1,\infty} \\ z_i =\left \{\begin{array}{ll} 1, p=0.5,\\ 0, p=0.5 \end{array}$$

Очевидно, реализовать такую дробь в компьютере невозможно, так как разрядная сетка компьютера ограничена. В компьютере можно формировать дискретные последовательности случайных чисел, которые не могут отличаться друг от друга только на величину меньше $$2^{-n}$$ ( $$n$$ - число разрядов в сетке компьютера). То есть непрерывного, "теоретического" распределения на компьютерах получить нельзя. Если эти числа равновероятны, то такое распределение случайных чисел называют квазиравномерным.

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

Случайная величина $$\xi $$, имеющая квазиравномерное распределение в интервале $$[0,1]$$, принимает значения

$$x_i=\cfrac{i}{2^n-1},\,\, i =0,1,2,\ldots,2^n-1$$

с вероятностями $$p_{i} =0.5^n$$.

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

$$M[\xi]=\cfrac{1}{2},\,\,D[\xi]=\cfrac{1\cdot(2^n+1)} {12\cdot(2^n-1)}.$$

Современные компьютеры имеют разрядность не менее 32. Следовательно, $$М[ \xi ] =М[ \gamma ]$$, а дисперсии тоже практически совпадают. Учитывая это, в дальнейшем квазиравномерное распределение будем называть равномерным и обозначать: $$\gamma \sim Rav[0, 1]$$.

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

аппаратный (физический);

табличный (файловый);

алгоритмический (программный).

Аппаратный способ. При этом способе случайные числа формируются специальным устройством. Источником случайных чисел чаще всего являются шумы в электронных приборах. Временные расстояния между шумовыми всплесками, превышающими подобранный уровень ограничения, фиксируются как случайные числа из распределения $$\gamma \sim Rav[0, 1]$$.

Преимущества такого способа:

количество случайных чисел неограниченно;

не требует затрат оперативной памяти;

требует малые вычислительные ресурсы компьютера.

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

трудность настройки;

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

обеспечение стабильности условий работы устройства - питания, влажности, температуры, старения приборов и элементов;

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

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

числа требуют однократную проверку при формировании или недоверии источнику;

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

Недостатки же очень существенны:

количество случайных чисел ограничено;

файл занимает место в оперативной памяти компьютера;

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

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

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

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

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

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

Недостатки:

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

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

В настоящее время практически везде применяются алгоритмические датчики случайных чисел (ДСЧ). Создание высокопроизводительных компьютеров существенно снижает роль первого недостатка (затраты машинного времени). Второй недостаток устраняется использованием в одной модели нескольких ДСЧ.

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

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

Качество алгоритмического датчика оценивается тем, насколько полно он удовлетворяет следующим требованиям:

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

числа должны быть статистически независимыми;

числа не должны повторяться;

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

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

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

подбирается начальное число $$x_{0}$$, например, $$x_{0} = 0.2152$$ ;

вычисляется квадрат $$x_{0}$$: $$x _{0}^{2} = 0,04631104$$ ;

первое случайное число $$x _{1} = 0.6311$$, то есть середина $$x _{0}^{2}$$ ;

вычисляется квадрат $$x _{1}$$: $$x _{1}^{2} = 0.39828721$$ ;

второе случайное число $$x_{2} = 0.8287$$ и т. д.

Такие ДСЧ теперь не используются: между числами имеется сильная корреляция, случайность отсутствует, при неудачно выбранном $$x_{0}$$ последовательность может быстро выродится, то есть $$x_{i} =0$$ при $$i \succ i^*$$.

В настоящее время очень широкое распространение в практике моделирования получил мультипликативный метод формирования случайной последовательности:

$$X_{i} _{+1} = \lambda X_i (modM),\,\, i = 0,1,\ldots ,$$

где $$X_{0} $$ - произвольное нечетное число, неотрицательное;

$$\lambda $$ - коэффициент, $$\lambda = 8t \pm 3$$, $$t $$ - любое целое положительное число;

$$M $$ - значение модуля. Для реализации на компьютере удобно $$M = p^q$$, где $$p $$ - основание системы счисления (2 или 10), $$q $$ - число разрядов в случайном числе.

В этом случае взятие числа по модулю сводится к выделению $$q $$ младших разрядов произведения $$\lambda X_{i}$$.

Алгоритм мультипликативного метода

Выбрать $$X_{0}$$, например, $$X_{0} = 234567.$$

Вычислить коэффициент $$\lambda $$. Пусть $$t = 1$$, тогда $$\lambda = 8?1-3 = 5$$.

Выбрать модуль $$M $$. Пусть система счисления десятичная ( $$p = 10$$ ), разрядность случайных чисел $$q = 6 . Тогда M = 10^6$$.

Вычислить произведение $$\lambda X_{0}$$: $$\lambda X_{0} = 5\cdot 234567 = 1172835 $$.

Найти остаток от деления по модулю $$M $$: $$X_{1} = \lambda X_{0}\cdot (mod10^6) = 1172835(mod10^6) = 172835 $$.

Найти число последовательности случайных чисел из интервала $$[0,1]$$: $$x_{1}=X_{1}/10^6 = 0.172835 $$.

Присвоить $$X_{0} = X_{1} $$ и перейти к п. 4.

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

Применяется и немного более сложный алгоритм:

$$X_{i} _{+1} = ( \lambda X_{i} +\mu )(modM), i = 0,1,\ldots ,$$

где $$\mu $$ - неотрицательное целое число.

Такой метод называется конгруэнтно-мультипликативным. При удачном подборе дополнительного параметра $$\mu $$ корреляция

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

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

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

В этом методе используется следующая теорема.

Теорема. Если случайная величина $$Y $$ имеет плотность распределения вероятностей $$f (y)$$, то распределение случайной величины

$$F(y)=\int\limits_{-\infty}^{y}{f(y)dy}$$

равномерно в интервале $$[0,1]$$, т. е.

$$F (y) = \gamma \sim Rav [0, 1] .$$

По определению, $$F(y) $$ является функцией распределения случайной величины $$Y $$.

Теорема может быть проиллюстрирована графиками, представленными на рис. 3.10.

Обозначим: $$x_{i} $$ - $$i$$ -е число из $$\gamma\sim Rav[0, 1]$$, $$y_{i} $$ - $$i$$ -е случайное число из произвольного распределения.

(рис 3.10) Иллюстрация к методу обратной функции Из (3.1) следует:$$x_{i} = \int\limits_{-\infty}{y_i}{f(y)dy} .$$

Моделировать равномерно распределенное случайное число $$x_{i }$$ мы уже умеем. Нужно найти неизвестное $$y_{i}$$, находящееся в верхнем пределе интегрирования.

Относительно $$y_{i} в$$ ыражение принимает вид:

$$y_{i}=F^{-1} ( x_{i}).$$

Отсюда и название - "метод обратной функции".

Пример 3.6. Получить формулу для моделирования случайных

чисел, распределенных по экспоненциальному закону, с параметром $$\lambda $$ (матожиданием $$\cfrac{1}{\lambda}$$ ).

Плотность $$f(y) $$ и функция $$F(y) $$ этого распределения имеют вид (рис. 3.11):

Решение

$$\xi =\int\limits_{-\infty}^{y_i}{f(y)dy}= \int\limits_{-\infty}^{0}{f(y)dy}+\int\limits_{0}^{y_i}{f(y)dy}= \int\limits_{0}^{y_i}{f(y)dy}=1-e^{-\lambda y_i}\\ x_i \in Rav[0,1];\,\, e^{-\lambda y_i}=1-x_i;\,\, {-\lambda y_i}=\ln{(1-x_i)};\,\, y_i=-\cfrac{1}{\lambda}{(1-x_i)}$$

Поскольку случайная величина $$(1-x_{i})$$ имеет равномерное распределение в интервале $$[0,1]$$, как и $$x_{i}$$, то справедливо:

$$y_{i} =\cfrac{1}{\lambda}\ln{x_{i}}.$$

Примеров подобного аналитического преобразования случайного числа $$x_{i} \in Rav[0, 1] $$ в случайное число из произвольного распределения немного, так как для многих законов распределения, встречающихся в практике моделирования, интеграл (3.1) относится к неберущимся, а численные методы решения увеличивают затраты машинного времени.

(рис 3.11) Плотность и функция экспоненциального распределения

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

Суть метода заключается в следующем.

Требуемый закон распределения случайной величины размещается в памяти компьютера в виде координат функции распределения. Каждая координата состоит из случайного числа $$y_{i} $$ и соответствующего значения функции распределения $$F ( y_{i} ):$$

$$y_{1}, F ( y_{1} ) / y_{2}, F ( y_{2} ) / \ldots / y_{k}, F ( y_{k} ) / \ldots / y_{n}, F ( y_{n} )$$

Чем больше координат, тем точнее будет моделирование. Приемлемая точность обеспечивается заданием 20…30 координат.

При обращении за очередным случайным числом нужного закона распределения сначала генерируется случайное число из $$Rav[0,1].$$ Это число сравнивается со значениями $$F ( y_{k} ), k =\overline{1, n}$$.

При совпадении выдается соответствующее случайное число $$y_{k} $$.

Если нет совпадения, то случайное число $$y_{k} $$ вычисляется из подобия треугольников, как показано на рис. 3.12.

(рис 3.12) Иллюстрация к методу кусочно-линейной аппроксимации

Из подобия треугольников ABC и AB'C' следует:

$$\cfrac{\Delta F(y_k)}{\Delta y_k} = \cfrac{x_i - F(y_k)}{y_i - y_k}$$

Отсюда по $$x_{i} \in Rav[0, 1]$$ находится значение $$y_{i}$$.

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

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

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

Теорема. Случайная величина

$$\eta=\cfrac{\sum\limits_{i=1}^{N}{x_i}-NM[x]}{\sqrt{ND[x]}}$$

где $$\sum\limits_{i=1}^{N}{x_i}$$ - сумма $$N $$ случайных чисел одного и того же распределения с матожиданием $$M[x]$$ и дисперсией $$D[x]$$ при $$N \to\infty $$ асимптотически стремится к нормальному распределению с $$M[\eta] = 0$$ и дисперсией $$D[\eta] = 1$$.

Удобно случайные числа $$x_{i} $$ брать из рассмотренного датчика $$у \sim Rav[0, 1]$$. В этом случае $$M[\gamma] = \cfrac{1}{2},$$ $$D[\gamma] =\cfrac{1}{12}$$.

Хорошее приближение к нормальному распределению получается уже при числе $$N = 6.$$ Каждое случайное число при $$N = 6$$ генерируется так:

$$y_i=\cfrac{\sum\limits_{i=1}^{N}{x_i-3}}{\sqrt{0.5}}=\sqrt{2}\left ( \sum\limits_{i=1}^{6}{x_i-3} \right )$$

Недостаток способа состоит в том, что он не экономичен, так как для генерирования одного случайного числа $$y_{i} $$ требуется шесть случайных чисел из распределения $$у \sim Rav[0, 1]$$.

В ряде случаев применяют датчики с числом $$N = 12$$. Тогда

$$y_i=\cfrac{\sum\limits_{i=1}^{12}{x_i-6}}{\sqrt{1}}=\sum\limits_{i=1}^{12}{x_i}-6 .$$

Если датчик случайных чисел нормального распределения выдает стандартную последовательность чисел с $$M = 0$$, $$\sigma = 1$$, то пересчет на произвольное значение характеристик выполняется так:

$$y_j =y_j\sigma + m,$$

где $$m $$ - требуемое значение матожидания;

$$\sigma$$ - требуемое значение среднего квадратического отклонения;

$$y^{'}_j$$ - случайное число из нормального распределения с математическим ожиданием $$m $$ и средним квадратическим отклонением $$\sigma$$.

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

3.6. Моделирование единичного события

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

Переход из одного состояния в другое - случайный. В любой момент времени система находится в одном состоянии с вероятностью $$P$$, в другом - с вероятностью $$1 - P $$.

Цель моделирования: имитировать состояние такого элемента.

Теорема. Пусть некоторое событие $$A $$ свершается с вероятностью $$p(A)$$. Это может быть отказ техники, поступление сообщения, уничтожение цели и т. п.

Моделью свершения такого единичного события A является попадание значения $$x_{i} $$ случайной величины $$\gamma$$, равномерно распределенной в интервале $$[0,1]$$, в числовой интервал $$[0, P(A)]$$.

Доказател ьство Как известно

$$P(a\le x \le b)=\int\limits_{a}^{b}{f(x)dx}$$

Для $$x_i\in \gamma \sim Rav[0,1], \, a=0,\, b=P(A)$$

$$P(0\le x_i \le P(A))=\int\limits_{0}^{P(A)}{f(x)dx}= \int\limits_{0}^{P(A)}{1\cdot dx} = P(A)$$

так как $$f(x)=1$$ для $$\gamma \sim Rav[0,1]$$ на интервале $$[0,P(A)]$$.

Пример 3.7. Пусть вероятность состояния элемента $$P(A) = 0.9.$$ В $$i $$ -ой реализации случайное число $$x_{i}\in\gamma\sim Rav[0, 1]$$ равно $$0.955.$$ Это означает, что в данной $$i$$ -ой реализации модели событие A не свершилось (рис. 3.13). Естественно, одна реализация ни о чем не говорит. Реальная ситуация будет отображена на множестве реализаций и чем их больше, тем точнее.

(рис 3.13) Событие A не произошло

Фрагмент алгоритма имитации в модели единичного случайного события приведен на рис. 3.14.

(рис 3.14) Фрагмент алгоритма имитации единичного события

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

Элемент системы (или система в целом) может находиться во многих (больше двух) несовместных состояниях. Известны вероятности нахождения системы в этих состояниях. Например, средство вооружения может находиться:

в боеготовом состоянии с вероятностью $$P_{1}$$ ;

в неисправном состоянии и ремонтироваться силами своего расчета, вероятность этого состояния $$P_{2}$$ ;

ремонтироваться в мастерской части - $$P_{3}$$ ;

ремонтироваться на заводе - $$P_{4}.$$ Очевидно, что $$P_{1} + P_{2} + P_{3} + P_{4} =1$$.

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

Алгоритм моделирования основан на следующей теореме.

Теорема. В полной группе несовместных событий моделью свершения события $$A_{m}$$, происходящего с вероятностью $$Р_{m}$$, является попадание значения $$х_i\in\gamma \sim Rav[0, 1]$$ в отрезок, равный $$Р_{m}$$, числовой шкалы $$\sum\limits_{m=1}^{n} Р_{m} = 1$$, где $$n$$ - число несовместных событий (рис. 3.15):

(рис 3.15) Событие A_{m} произошло

Доказател ьство

Введем численные обозначения концов отрезков $$Р_{m}$$ по нарастанию:

$$l_r = \sum\limits_{m=1}^{r} Р_{m}$$

В этом случае, согласно теореме, условием свершения события $$A_{m} $$ является:

$$l_{m-1} \prec x_i \le l_m.$$

Следовательно

$$P(l_{m-1} \prec x_i \le l_m) = \int\limits_{l_{m-1}}^{l_m}{\gamma(x)dx} = \int\limits_{l_{m-1}}^{l_m}{1\cdot dx} = \left ( l_m = \sum\limits_{1}^{m}{P_m}\right ) - \left ( l_{m-1} = \sum\limits_{1}^{m-1}{P_m}\right ) = P_m$$

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

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

Первый вариант (рис. 3.16а) применяется тогда, когда число возможных исходов невелико и не равно степени по основанию два.

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

(рис 3.16) Варианты алгоритма определения исходов по жребию

Третий вариант алгоритма (рис. 3.16в) в цикле определяет исход (событие), номер которого присваивается переменной $$k $$. Далее этот номер используется для организации нужной работы алгоритма. Применение данного алгоритма будет показано в главе 6 (п. 6.7 и п. 6.8).

Пример 3.8. Канал передачи данных может находиться в одном из четырех несовместных состояниях:

$$A_{1} $$ - исправен и свободен, $$P_{1} = 0,15$$ ;

$$A_{2}$$ - исправен и занят, $$P_{2} = 0,4$$ ;

$$A_{3}$$ - неисправен, $$P_{3} = 0,25$$ ;

$$A_{4}$$ - подавлен помехами, $$P_{4} = 0,2.$$

Решение

Представим необходимые для определения исходов по жребию данные табл. 3.3.

Данные для определения исходов по жребию

Вероятности	Событие
	$$A_{1}$$	$$A_{2}$$	$$A_{3}$$	$$A_{4}$$
Вероятности событий	0,15	0,4	0,25	0,2
Суммарные вероятности ( $$l_{r}$$ )	0,15	0,55	0,8	1,0
Номера интервалов ( $$r$$ )	1	2	3	4

Предположим, что при выполнении $$i$$ -ой реализации датчик равномерно распределенных случайных чисел $$\gamma\sim Rav[0, 1]$$ сгенерировал $$x_{i} = 0.525.$$ Путем последовательных сравнений определяется, что $$l _{1} < 0,525 \le l_{2}$$. Значит в данной реализации канал находится в состоянии $$A _{2}$$ - исправен и занят.

3.8. Моделирование совместных независимых событий

Рассмотрим моделирование совместных независимых событий.

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

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

Пусть независимые события A и B происходят с вероятностями $$P(A) $$ и $$P(B)$$ соответственно. Например, это могут быть отказы монитора и процессора компьютера.

Моделирование такой ситуации может быть выполнено двумя способами:

определение совместных исходов выбором по жребию;

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

3.8.1. Определение совместных исходов по жребию

Прежде всего, по вероятностям $$P(A)$$ и $$P(B) $$ нужно определить вероятности возможных исходов, т. е. появления совместных независимых событий. Возможные исходы совместного события $$Q_{i} $$ и соответствующие вероятности $$P_{i} $$ представлены в табл. 3.4.

Возможные исходы совместного события

$$Q_i$$	$$AB$$	$$\overline{A}B$$	$$A \overline{B}$$	$$\overline{A}\overline{B}$$
$$P_i$$	$$P(A)P(B)$$	$$[1-P(A)]P(B)$$	$$P(A) [1-P(B)]$$	$$1-P_1$$
$$l_r$$	$$l_1= P(A)P(B)$$	$$l_2= l_1+[1-P(A)]P(B)$$	$$l_3=l_2+ P(A) [1-P(B)]$$	$$1$$

Совместное событие в $$i $$ -ой реализации определяется выбором исхода по жребию.

Если случайное число $$x_{i} \in \gamma \sim Rav[0, 1]$$ при очередной реализации окажется, например, на участке $$l _{1} \prec x_{i} \le l_{2}$$, то в данной реализации фиксируется свершение сложного события $$AB $$. Если же окажется $$x_{i} > l_{3}$$, то фиксируется событие $$AB $$. Алгоритм может быть построен по одному из приведенных на рис. 3.16 вариантов.

3.8.2. Последовательная проверка исходов

Алгоритм способа последовательной проверки исходов приведен на рис. 3.17.

(рис 3.17) Алгоритм последовательной проверки исходов

Проверку свершения каждого из совместных событий надо осуществлять разными случайными числами, так как события независимые. При первом способе достаточно одного случайного числа $$x_{i}$$, но сравнений может быть больше. Кроме того, нужно предварительно рассчитывать вероятности возможных исходов.

3.9. Моделирование совместных зависимых событий

Пусть события A и B имеют вероятности свершения $$P(A) $$ и $$P(B)$$ соответственно. Условная вероятность $$P(B/ A) $$ известна.

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

Пример 3.9. При испытании нового автомата определены вероятности горизонтального и вертикального отклонений пробоин от точки прицеливания $$P(A) = P(r _{1} \le 10см )$$ и $$P(B) = P(r_{2}\le 10см ).$$

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

$$P(B/\overline{A})=\cfrac{P(B) - P(A)P(B/A)}{1-P(A)}$$

Соответствующий фрагмент модели приведен на рис. 3.18.

(рис 3.18) Алгоритм моделирования совместных зависимых событий

Пример 3.10. В ремонтное подразделение поступают вышедшие из строя средства связи (СС). В каждом СС могут быть неисправными в любом сочетании блоки A, B, C . Вероятности выхода из строя блоков $$P_{A}$$, $$P_{B}$$, $$P_{C} $$ соответственно. Ремонт производится путем замены неисправных блоков исправными блоками. В момент поступления неисправного СС вероятности наличия исправных блоков $$P_{HA}$$, $$P_{HB} $$, $$P_{HC} $$ соответственно. При отсутствии хотя бы одного из исправных блоков A, B, C ремонт неисправного СС не производится.

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

Решение

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

Вероятности появления неисправных блоков

$$Q_i$$	$$ABC$$	$$AB\overline{C}$$	С другими блоками
$$P_i$$	$$P_{A}P_{B}P_{C}$$	$$P_{A} P_{B} (1- P_{C} )$$
$$l_r$$	$$l_{1} =P_{A} P_{B} P_{C}$$	$$l_{2} = l_{1} + P_{A} P_{B} (1- P_{C} )$$	1

Вероятности наличия исправных блоков

$$Q_{Hi}$$	$$ABC$$	$$AB\overline{C}$$	С другими блоками
$$P_{Hi}$$	$$P_{HA} P_{HB} P_{HC}$$	$$P_{HA} P_{HB} (1- P_{HC})$$
$$l_{Hr}$$	$$l_1 = P_{HA} P_{HB} P_{HC}$$	$$l_{2} =l_{1} + P_{HA} P_{HB} (1- P_{HC})$$	1

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

Алгоритм имитационной модели приведен на рис. 3.19.

В алгоритме приняты следующие обозначения:

$$N_{0} $$ - заданное количество реализаций модели;

$$N $$ - счетчик количества реализаций модели; $$M $$ - счетчик числа отремонтированных СС за $$N $$ реализаций модели;

$$A$$ - абсолютное количество отремонтированных СС;

$$D $$ - относительное количество отремонтированных СС.

(рис 3.19) Алгоритм модели функционирования системы ремонта

Согласно постановке задачи в блоках 3…7 по данным табл. 3.5 разыгрывается, с какими неисправными блоками поступает СС в ремонт. В результате розыгрыша определяется номер интервала (столбца табл. 3.5) и запоминается в переменной $$k $$.

Аналогично в блоках 8…11 разыгрывается по данным табл. 3.6 наличие в ремонтном подразделении необходимых блоков для замены.

Если такие блоки имеются, т. е. выполняется условие $$k =i $$ в блоке 12, в счетчик $$M $$ (блок 13) добавляется единица.

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

Случайная величина $$X (t ),$$ зависящая от одного неслучайного вещественного аргумента $$t $$, называется случайным процессом. $$X (t )$$ является случайной величиной при каждом фиксированном

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

Определим два понятия, присущие случайным процессам: сечение и реализация (рис. 3.20).

Сечением случайного процесса $$X (t )$$ называется случайная величина $$x(t_{j} ),$$ являющаяся значением случайного процесса в фиксированный момент времени $$t_{j} $$.

Реализацией случайного процесса $$X (t )$$ называется функция времени $$x_{i} (t ),$$ описывающая течение процесса в некотором $$i $$ -м опыте.

(рис 3.20) Реализации и сечения случайного процесса

Случайный процесс $$X (t )$$ и аргумент $$t $$ могут быть дискретными или непрерывными.

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

разбить интервал исследования на M временных точек $$t_j (j = \overline{1, M})$$, которых должно быть столько, чтобы обеспечить необходимую точность воспроизведения исследуемого процесса;

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

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

Случайные процессы могут быть:

стационарные;

нестационарные.

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

$$M[x(t)]=M[x];\,\, D[x(t)] = D[x].$$

Такие процессы называют стационарными.

Что касается автокорреляционной функции $$K_{k}(t_{k},t_{3})$$, то ее значение в стационарном процессе зависит только от разности $$t_{k}-t_{3} $$ и не зависит от того, в каком месте временной оси находятся точки $$t_{k} $$ и $$t_{3}$$.

Для стационарного процесса нет необходимости определять искомые характеристики для всех $$M $$ сечений, а достаточно только для одного сечения $$N $$ реализаций случайного процесса. То есть вместо $$M\times N $$ измерений выполнить только $$N $$ измерений. По данным этих измерений рассчитываются оценки $$x(t_j) $$ и $$S^2_{x(t_j )}$$, которые в силу стационарности и являются оценками характеристик всего случайного процесса $$M[x]$$ и $$D[x]$$.

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

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

Важнейшим свойством случайного процесса является свойство эргодичности.

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

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

$$\cfrac{1}{N}\sum\limits_{i=1}^{N}{x_i(t_j)}= \cfrac{1}{M}\sum\limits_{j=1}^{M}{x_i(t_j)}$$

Свойством эргодичностиобладают многие случайные процессы и, в том числе, все стационарные.

Таким образом, можно сформулировать определение эргодиче-ского процесса.

Случайный процесс $$X (t )$$ называется эргодическим, если его основные характеристики $$M [x]$$ и $$D[x]$$ могут быть получены не только усреднением по множеству реализаций, но и усреднением по времени одной реализации.

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

3.11. Способы продвижения модельного времени

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

реальное время системы, функционирование которой имитируется;

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

машинное время имитации, отражающее затраты ресурса времени компьютера.

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

Следовательно, время в модели, то есть модельное время (МВ), продвигается дискретно, скачками.

Продвижение времени в модели может быть организовано двумя способами:

продвижение модельного времени с фиксированным переменным шагом $$\Delta t $$ ;

продвижение модельного времени до очередного события (по принципу $$\Delta x $$ ).

Сущность первого способа поясним временными диаграммами, показанными на рис. 3.21.

На диаграммах а…г показаны моменты смены дискретных состояний элементами 1…4 системы. На диаграмме д - временная последовательность смены состояний системой. На диаграмме е - точки модельного времени, то есть время смены состояний системы, показанных на диаграмме д.

Так как моменты модельного времени на диаграмме е не связаны с моментами появления событий а…г, то имитационная модель с фиксированным шагом продвижения времени искажает действительные процессы в системе: разновременные события представляются одновременными, моменты свершения событий фиксируются, как правило, с опозданием. Уменьшая величину $$\Delta t $$, можно уменьшить искажение действительного процесса. Однако это приводит к увеличению затрат машинного времени, особенно, если интервалы между сменами состояний в среднем больше, чем $$\Delta t $$.

На диаграмме ж рис. 3.20 демонстрируется сущность второго способа. Она заключается в том, что модельное время сдвигается вперед не на фиксированную величину $$\Delta t $$, а точно до времени наступления самого раннего из очередных событий - на $$\Delta x$$.

Видно, что недостатки, присущие первому способу, здесь исключены: события рассматриваются и моделируются в моменты их свершения, и одновременно (события a₁,г₁), если у них одинаковое время появления. Промежутки времени, когда в модели "ничего не происходит", пропускаются без особых затрат машинного времени. Эти пропуски все равно учитываются в модельном времени.

(рис 3.21) Временная диаграмма работы модели

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

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

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

Однако встречаются имитационные модели, которые предназначены для работы в реальном (естественном) масштабе времени. Это, например, некоторые типы тренажеров, работающие в интерактивном режиме с человеком. Заметим, что и в этом случае модельное время продвигается скачками - способом $$\Delta t $$ или $$\Delta x $$. Следовательно, возникает проблема синхронизации модельного времени с естественным временем.

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

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

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

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

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

Пример 3.11. На узел связи поступают заявки на передачу сообщений. Интервалы времени поступления заявок подчинены показательному закону с математическим ожиданием $$Т_{1} $$. На узле связи имеются два канала передачи данных. При поступлении очередной заявки в интервале времени $$[0\ldots T_{2}]$$ вероятности того,

что каналы А и В будут свободны, соответственно равны $$Р_{1А}$$ и $$Р_{1В}$$. При поступлении заявок после времени $$Т_{2} $$ вероятности того, что каналы А и В будут свободны, соответственно равны $$Р _{2 А } $$ и $$Р _{2В}$$. Сообщение передаётся по любому свободному каналу. Если оба канала заняты, заявка теряется.

Построить алгоритм имитационной модели "Обработка запросов на узле связи" с целью определения абсолютного и относительного числа обслуженных заявок из их общего количества, поступивших на узел связи за время $$T_{MOD}$$, $$T_{MOD} \succ Т_{2}$$.

Решение

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

Алгоритм модели приведен на рис. 3.22.

(рис 3.22) Алгоритм модели "Обработка запросов на узле связи"

В алгоритме модели:

$$S $$ - суммарное количество заявок, поступивших за $$n $$ реализаций модели;

$$t $$ - текущее модельное время;

$$N $$ - заданное количество реализаций модели;

$$n$$ - счетчик текущего числа реализаций модели;

$$М $$ - суммарное число обслуженных заявок за $$n$$ реализаций;

$$A $$ - абсолютное число обслуженных заявок;

$$Y $$ - относительное число обслуженных заявок.

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

внутренний определяется временем моделирования $$T_{MOD}$$ (блок 15), а второй, внешний - количеством реализаций модели $$N $$ (блок 17).

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

Далее мы рассмотрим способы, позволяющие в некотором роде унифицировать построение моделей сложных систем с продвижением в них модельного времени по принципу $$\Delta t $$ или $$\Delta х $$, а также до ближайшего события.

3.12. Модель противоборства двух сторон

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

Для начала, чтобы не загромождать алгоритм модели, введем два ограничения:

вероятность поражения цели одним выстрелом $$\le 0.1$$ ;

вероятность одновременного поражения двумя средствами друг друга пренебрежимо мала.

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

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

Учебная задача состоит в том, чтобы изучить структуру конкретного алгоритма, в котором продвижение модельного времени реализовано с фиксированным шагом $$\Delta t $$.

Введем обозначения:

$$А, В $$ - первоначальные численности группировок;

$$a, b$$ - текущие значения численностей группировок $$А $$ и $$В $$ соответственно;

$$a_{k}, b_{k} $$ - число оставшихся средств каждой из сторон в конце $$k$$ -го интервала моделирования (после $$k$$ -ой реализации модели);

$$n $$ - текущее число реализаций модели;

$$N $$ - заданное число реализаций модели случайного процесса;

$$t $$ - текущее время;

$$T $$ - длительность интервала моделирования;

$$\overline{a}(T ), \overline{b}(T)$$ - средние за $$N $$ реализаций модели численности оставшихся средств сторон в конце каждого интервала моделирования:

$$\overline{a}(T)=\cfrac{\sum\limits_{k=1}^{N}{a_k}}{N},\,\, \overline{b}(T)= \cfrac{\sum\limits_{k=1}^{N}{b_k}}{N},$$

$$i, j $$ - переменные счета средств сторон $$А $$ и $$В $$ соответственно .

Идея построения имитационной модели состоит в реализации модели противоборства $$N $$ раз и фиксации остатков сторон после каждой реализации. По выполнении $$N $$ реализаций будут определены оценки численностей сторон $$a(T), b(T)$$.

Определим величину $$\Delta t $$. На выбор величины $$\Delta t $$ влияют два противоречивых требования:

на протяжении отрезка времени $$\Delta t $$ не должно происходить много событий, так как они будут зафиксированы как одновременные, что исказит исследуемый реальный процесс;

на протяжении отрезка времени $$\Delta t $$ должно произойти хотя бы одно событие, иначе будет много "пустых" прогонов модели, что увеличит машинное время.

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

Например, скорострельность средств поражения стороны $$A$$

$$\lambda _{A} = 2выстр./мин $$, а средств стороны $$B- \lambda _{B} = 3выстр./мин $$. Так как $$\lambda_{B}\succ \lambda_{A}$$, то $$\Delta t$$ следует выбрать исходя из условия $$\Delta t \prec 1 / \lambda_{B} $$.

Алгоритм имитации противоборства сторон состоит из четырех модулей:

$$М1$$ - установка начальных условий;

$$М2 $$ - продвижение модельного времени;

$$М3$$ - формирование результата заданной точности;

$$М4$$ - имитация противоборства сторон.

Блок-схема алгоритма имитации противоборства двух сторон показана на рис. 3.23.

Модуль установки начальных условий ( $$М1$$ ). Он состоит из трех блоков 1… 3.

Блок 1 - установка начальных условий на весь процесс моделирования: число реализаций $$N $$ модели, интервал исследования $$T$$, величина временного шага $$\Delta t$$, установка в нуль ячеек $$\sum{a_k},\sum{b_k}$$

Блок 2 - установка начальных условий на очередную реализацию процесса: восстановление численностей сторон $$a= A,b = B $$ и исходного времени $$t = 0 $$.

Блок 3 - установка начальных условий на очередной отрезок $$\Delta t $$ модельного времени: подготовка перебора $$(i = 0, j = 0)$$ средств

поражения каждой стороны.

Модуль продвижения модельного времени ( $$М2$$ ). Состоит из двух блоков 4..5.

Блок 4 - продвижение модельного времени на очередной временной отрезок $$\Delta t (t = t + \Delta t)$$.

Блок 5 - проверка условия окончания очередной реализации модели $$(t \ge T) $$.

Модуль формирования результата и обеспечения заданной точности ( $$М3$$ ). Состоит из блоков 6…9.

Блок 6 - накопление суммы остатков средств каждой стороны за текущее количество интервалов моделирования (реализаций модели).

(рис 3.23) Блок-схема алгоритма модели противоборства двух сторон

Блок 7 - счетчик числа реализаций модели ( $$n = n+1$$ ) .

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

Блок 9 - формирование конечного результата моделирования.

В данном случае - $$\overline{a} $$ и $$\overline{b}$$.

При необходимости определяются и другие статистические характеристики.

Модуль имитации противоборства сторон ( $$М4$$ ). Этот модуль - основной. Если структура предыдущих модулей, в общем-то, стандартна, то реализация данного модуля носит функционально-индивидуальный характер. Модуль включает блоки 10…25. Блок-схема алгоритма модуля представлена на рис. 3.24.

Каждое средство противоборствующих сторон идентифицируется его номером. Номера средств стороны $$A$$ - $$i = \overline{1,A} $$, средств стороны $$B $$ - $$j = \overline{1,B}$$.

Блок 10 - выбор очередного средства стороны $$A$$: $$i = i + 1$$.

Блок 11 - проверка: все ли средства стороны $$A $$ получили право на выстрел? Если $$i \succ A $$, то управление передается блоку 18 для имитации выстрелов средствами стороны $$B$$. В противном случае управление передается блоку 12.

Блок 12 - проверка: боеспособно ли выбранное средство? Состояние средств сторон $$A$$ и $$B $$ определяют переменные $$\alpha _{i} $$ и $$\beta _{j} $$:

$$\left\{ \begin{array}{l} \alpha_i =0, \text{ если } i\text{-ое средство стороны }A \text{ уничтожено;}\\ \alpha_i =1, \text{ если } i\text{-ое средство стороны }A \text{ боеспособно;} \end{array} \\$$ $$\left\{ \begin{array}{l} \beta_j =0, \text{ если } j\text{-ое средство стороны }B \text{ уничтожено;}\\ \beta_j =1, \text{ если } j\text{-ое средство стороны }B \text{ боеспособно;} \end{array}$$

Если окажется $$\alpha _{i} = 0$$, то управление передается блоку 10 для

выбора очередного средства стороны $$A $$. Иначе - переход к блоку 13.

Блок 13 - выбор цели из средств стороны $$B $$. Выбор цели может быть организован либо случайной, либо детерминированной процедурами. Самый простой способ: последовательная проверка средств стороны $$B $$ с выбором первого непораженного средства.

Блок 14 - проверка выбранной цели: не уничтожена ли она была на предыдущих этапах данной реализации модели? Если $$\beta_j =0$$, то переход к блоку 13 для выбора непораженной цели.

Иначе - переход к блоку 15 для имитации выстрела.

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

(рис 3.24) Блок-схема алгоритма модуля противоборства двух сторон

Блок 15 - обращение к ДСЧ за равномерно распределенным случайным числом $$х_i \sim Rav[0, 1]$$.

Блок 15.1 - проверка результата выстрела. Если $$х_i \le Р_{аb}$$, цель поражена и управление передается блоку 16 для фиксации этого факта. Если $$х_{i} \succ Р_{ab}$$, то промах и управление передается блоку 10 для выбора очередного стреляющего средства стороны $$A $$.

Блок 16 - уменьшение числа средств стороны $$В$$: $$b=b-1 $$ и установка признака состояния пораженного средства: $$\beta_j = 0$$.

Блок 17 - не уничтожена ли вся группировка $$В $$? Если да, то данная реализация модели заканчивается и управление передается блоку 6 для фиксации оставшихся боеспособных средств стороны $$A$$. Если нет, то управление передается блоку 10 для выбора очередного средства стороны $$A $$ и т. д.

После предоставления права на выстрел всем средствам стороны $$A $$ соответствующее право дается средствам стороны $$В $$ - переход из блока 11 в блок 18.

Функции блоков 18… 25 попарно одинаковы с функциями блоков 10…17, изменены только обозначения - вместо $$A $$ указано $$В$$, вместо $$\alpha _{i} \to \beta_j $$, вместо $$Р_{ab} \to Р_{ba}$$ и т. д.

По окончании перебора всех средств стороны $$В $$ и, если не зафиксировано полное уничтожение средств стороны $$A (a = 0)$$, управление передается блоку 3 для моделирования очередного скачка времени на величину $$\Delta t$$.

Примерная диаграмма изменения численностей сторон $$A$$ и $$В $$ в $$k$$ -й реализации на интервале $$0\ldots Т $$ показана на рис. 3.25.

(рис 3.25) Диаграмма изменения численности сторон

3.13. Модель противоборства как процесс блуждания по решетке

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

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

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

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

Противоборствуют две стороны: $$A$$ и $$B$$.

Первоначальные численности сторон: $$A $$ и $$B $$.

Текущие численности сторон: $$a $$ и $$b $$ ; $$a = \overline{0, A}$$, $$b = \overline{0, B}$$.

Время "жизни" каждого средства случайно, имеет экспоненциальное распределение с параметрами $$\lambda_{a}$$ и $$\mu_{b}$$ сторон $$A$$ и $$B$$.

Обе стороны образуют единую систему $$S_{ab} $$. В ходе боя численности сторон изменяются (уменьшаются) или остаются неизменными. Мы это трактуем как переход системы из одного состояния в другое: $$S_{a,b} \to S_{a-1,b}$$ или $$S_{a,b} \to S_{a,b-1}$$ и т. д.

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

Граф состояний приведен на рис. 3.26.

(рис 3.26) Граф состояний моделируемой системы противоборства Исходное состояние системы S_{a,b} .

Стороны ведут взаимоуничтожающий огонь. В предположении ординарности потока событий, переводящего систему из одного состояния в другое, система из состояния $$S_{a,b}$$ может перейти только лишь в состояние $$S_{a-1,b}$$ или $$S_{a,b-1}$$ (или остаться в прежнем состоянии $$S_{a,b}$$ ). Очевидно, направления переходов в графе - решетке только вправо или вниз.

Процесс противоборства заканчивается при полном уничтожении средств одной из сторон, т. е. при переходе в состояние $$S_{0,b}$$ или $$S_{a,0},$$ а также по окончании временного интервала исследования $$T $$. Так заканчивается одна, $$k $$ -я реализация процесса противоборства. Зафиксировав остатки средств сторон $$a_{k} $$ и $$b_{k} $$ и усреднив их за $$N $$ реализаций, найдем решение поставленной задачи - определение оценок численностей сторон на момент окончания операции.

Продвижение времени в модели - фиксированными шагами $$\Delta t$$. Величина $$\Delta t $$ такова, чтобы на этом временном участке происходило не более одного события.

Суммарные интенсивности огня сторон $$A$$ и $$B $$ зависят от численностей их боеспособных средств и равны $$a\lambda_{a} $$ и.

Следовательно, вероятности уничтожения одного из средств (перехода системы в очередное состояние) равны:

$$P(S_{a,b}\to S_{a,b-1}) \cong a\lambda_a\Delta t = P_a, \\ P(S_{a,b}\to S_{a-1,b}) \cong b\mu_b\Delta t = P_b$$

Эти приближенные равенства тем точнее, чем меньше $$\Delta t$$.

Блок схема алгоритма ИМ противоборства двух сторон, построенная способом "блуждания частицы по решетке" представлена на рис. 3.27.

Блоки 1 и 2 - установка начальных условий на весь процесс моделирования и на каждую очередную реализацию.

Блок 3 - расчет вероятностей $$P_{a}$$ и $$P_{b}$$.

Блоки 4…10 - определение исхода противоборства на очередном временном отрезке $$\Delta t$$. Если $$x_{i} \le P_{a}$$, то свершилось событие $$S_{a,b-1}$$, иначе - проверка условия $$x_{i} \le P_{b}$$? Если это условие выполняется, то свершилось событие $$S_{a-1,b}$$, если не выполняется, то пораженных средств не оказалось. Вероятность этого события $$S_{a,b}\to S_{a,b}$$ равна $$1-P_{a}-P_{b}$$.

Блоки 11… 16 аналогичны соответствующим блокам предыдущей имитационной модели противоборства двух сторон (см. рис. 3.22).

Сравним обе рассмотренные модели противоборства двух сторон.

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

Однако, первую модель несложно развить, например, учесть неоднородность средств в каждой группировке, указав соответствующие значения $$P_{ab} $$ и $$P_{ba} $$.

Последняя модель таких возможностей не имеет.

(рис 3.27) Блок-схема модели "блуждание частицы по решетке"

3.14. Типовая схема имитационной модели с продвижением времени по событиям

Уточним понятие события.

Под событием будем понимать смену состояния системы.

Виды событий:

выход из строя технических средств;

восстановление работоспособности техники;

поступление сообщения на узел связи;

начало передачи сообщения;

конец передачи сообщения;

уничтожение цели и т. п.

Вид события связан с типом элемента, сменяющего свои состояния. Смена состояний системой есть объединение смен состояний ее элементов.

Будем считать, что событие совершается в конкретный момент времени мгновенно.

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

Каждому событию соответствует пространственно-временная точка $$t ^{j,k}$$, где $$t $$ - момент свершения события; $$k $$ - тип элемента системы, сменившего состояние; $$j $$ - вид или номер события этого элемента.

Различают события активные и пассивные.

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

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

Пассивное событие - событие, возникшее под воздействием активного события.

Например:

начало передачи сообщения по каналу связи, так как оно зависит от другого активного состояния - поступления сообщения на узел связи;

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

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

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

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

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

Обычно элемент, в котором происходит активное событие, называют активным, пассивное событие - пассивным.

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

Множество точек $$t^{j,k}$$ представляет собой поведение системы во времени. Множество точек $$t^{j,k}$$ активных событий называется списком событий.

Список событий может быть сформирован либо перед началом моделирования, либо формироваться в ходе моделирования. На рис. 3.20 список событий представлен диаграммой ж.

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

модуль установки начальных условий;

модуль продвижения модельного времени;

модули реакции;

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

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

Модуль установки начальных условий.

Модуль устанавливает значения разного рода констант, например, заданное число $$N $$ реализаций модели, критерий завершения интервала моделирования $$T $$ и др., а также элементы списка событий $$t^{j,k}$$ - все или только начальные для каждого вида события.

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

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

ближайшего события $$\min t^{j,k}$$, и фиксация очередного момента модельного времени:

$$t = t_{0} + \min t^{j,k}.$$

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

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

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

подсчет числа свершившихся событий данного вида;

прогноз времени свершения $$t^{j,k}$$, очередного события данного вида и занесение этого времени в список событий (если список событий формируется в ходе моделирования). Прогноз осуществляется с помощью соответствующего датчика случайных чисел, адекватно имитирующего временные интервалы между событиями.

Модуль обеспечения заданной точности и достоверности ведет подсчет числа реализаций модели $$n = n +1$$. При достижении $$n\succ N $$ моделирование заканчивается.

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

Взаимодействие модулей показано на рис. 3.28.

Переход на соответствующий модуль реакции $$A$$, $$B $$ и т. д. происходит по значению индекса $$k $$ текущего значения $$t^{j,k}$$.

(рис 3.28) Блок-схема модели с продвижением времени по событиям

Безусловный выход из любого модуля реакции - в блок выбора очередного значения $$\min t^{j,k}$$ модуля продвижения времени.

В общий цикл предлагаемой структуры моделирующего алгоритма удобно вписывается определение конца интервала исследования $$T $$ и момента конца моделирования ( $$n = N $$ ) . Для этого вся система интерпретируется как один обобщенный активный элемент, имеющий два состояния: включен и выключен.

В списке событий этому обобщенному элементу соответствует одна пространственно-временная точка $$t^{j,k}= t_{0} +T $$. Здесь $$t_{0}$$ - начальная точка отсчета интервала исследования, обычно $$t_{0} =0$$. Когда в процессе работы моделирующего алгоритма окажется, что $$\min t^{j,k}=t_{0}+T $$, управление будет передано специальному модулю реакции $$\Theta$$, который зафиксирует конец текущей реализации $$(n = n + 1)$$.

Этот модуль также в случае $$n\prec N $$ передаст управление блоку ( $$НУ(N)$$ ) на выполнение очередной реализации, а в случае выполнения $$n = N $$ - модулю формирования результатов.

Ранее было отмечено, что некоторые элементы в зависимости от состояния могут быть активными либо пассивными.

В случае перехода такого элемента в пассивное состояние надо

исключить возможность выбора его блоком "Выбор $$\min t^{j,k}$$ ". Для этого в список событий записывается время $$t^{j,k} \succ T $$ его переключения на свой модуль реакции. В этом случае модуль продвижения времени никогда не найдет этой временной точки, так как интервал исследования будет заканчиваться раньше.

Но когда этот элемент вновь станет активным, время его активного события $$t^{j,k}$$, должно быть занесено в список событий.

3.15. Имитационная модель системы массового обслуживания

Рассмотрим имитационную модель СМО следующего вида:

многофазная;

многоканальная;

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

разомкнутая;

абсолютно надежная;

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

СМО такого вида показана на рис. 3.29.

На рис. 3.29:

$$\Gamma_{1}\ldots \Gamma_{G}$$ - источники потоков заявок на обслуживание;

$$K _{1,Q} \ldots K_{S,Q_S}$$ - каналы обслуживания;

$$S $$ - число фаз;

$$Q_1\ldots Q_S $$ - число каналов в каждой фазе;

$$L_{1}\ldots L_{S} $$ - очереди заявок на входах соответствующих фаз.

(рис 3.29) Многофазная СМО

Такой СМО может быть, например, сборочный цех радиоэлектронного предприятия. Из нескольких цехов ( $$\Gamma_{1}\ldots \Gamma_{G}$$ ) поступают готовые блоки для сборки изделия. Сначала они поступают на диспетчерские посты, на которых происходит оформление соответствующей документации. После этого блоки направляются на сборочные участки.

Диспетчеризация - первая фаза обслуживания.

Сборка изделий - вторая фаза обслуживания.

Затем идет стендовый контроль и регулировка - это третья фаза обслуживания.

Четвертая фаза - приемка готового изделия (например, военная приемка)

Далее может быть упаковка изделия - пятая фаза.

Упакованные изделия поступают в транспортное подразделение для отправки - это шестая фаза.

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

Вернемся к общей постановке задачи.

Целью моделирования является определение показателей исхода операции массового обслуживания на временном интервале $$0\ldots T $$. Например, оценка относительной пропускной способности СМО, среднее количество заявок в очередях, коэффициент использования каналов в каждой фазе и др.

Проанализируем элементы СМО.

Входные потоки заявок на обслуживание

Поступление заявок на вход первой фазы - событие активное.

Следовательно, в модели должны быть предусмотрены модули реакции на это событие. Каждый поток обслуживает свой модуль реакции, значит, модулей реакции этого типа должно быть столько, сколько имеется потоков заявок, то есть $$G $$. Моменты поступления заявок имитируются соответствующими генераторами (датчиками) случайных интервалов между заявками потоков $$\Gamma_{1}\ldots \Gamma_{G}$$.

Каналы обслуживания

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

$$Q_{1} +Q_{2} +\ldots +Q_{S} .$$

Очереди заявок на обслуживание на входах фаз

Эти элементы СМО - пассивные, число состояний неограниченно.

Очереди свободных каналов в каждой фазе обслуживания

Эти элементы СМО также пассивные, число состояний равно числу каналов в каждой фазе.

Обобщенный элемент "СМО в целом"

Характеризуется одним активным событием - завершение интервала исследования $$T $$. Следовательно, для обработки этого события в модели предусматривается один специальный модуль реакции.

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

$$R = G + \sum\limits_{i=1}^{S}{Q_i} + 1$$

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

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

Таким образом, имитационная модель СМО будет иметь всего четыре модуля реакции:

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

модуль реакции на освобождение канала не последней фазы (В);

модуль реакции на освобождение канала последней фазы (С);

модуль реакции на завершение интервала исследования $$T $$ - реакция на изменение состояния обобщенным элементом "СМО в целом".

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

Блок 1 - установка начальных условий на весь процесс моделирования. Вводятся константы $$N$$, $$T $$, указываются число источников заявок, число каналов. Счетчик числа реализаций $$n $$ устанавливается в нуль (или в $$N $$ - в зависимости от того, как организован счет реализаций модели). Указываются данные планирования эксперимента с моделью и др.

Блок 2 - установка начальных условий на очередную реализацию модели. Текущее время, счетчики числа поступивших и обслуженных заявок устанавливаются в нуль. Каналы обслуживания и очереди к ним устанавливаются в исходные состояния (установки исходных состояний определяются требованиями к исследованию) и др.

(рис 3.30) Блок схема имитационной модели СМО

Наборы исходных данных $$НУ(М)$$ и $$НУ(n)$$ определяются конкретным назначением моделируемой СМО и характером последующих экспериментов.

Блоки 3, 4 реализуют продвижение времени по событиям.

Блок 5 - выбор модуля реакции. В зависимости от индекса $$k$$ пространственно временной точки $$\min t^{j,k}$$ передает управление соответствующему модулю реакции - $$A$$, $$B$$, $$C $$ или $$\Theta.$$

В рассматриваемой модели элементы СМО имеют только два состояния, одно из которых является начальным. Поэтому индекс $$j $$ не имеет смысла.

Рассмотрим состав и функционирование модулей реакции.

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

Блок-схема алгоритма модуля реакции $$А$$ приведена рис. 3.31.

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

(рис 3.31) Блок-схема алгоритма модуля реакции А

Блок А6 - проверка наличия свободных каналов первой фазы. Количество свободных каналов в очереди обозначено: $$l_{ск} $$.

Если все каналы заняты $$(l_{ск} = 0)$$, управление передается блоку А7. Если свободные каналы есть $$(1_{ск} \succ 0)$$, управление передается блоку А9.

Блок А7 - размещение поступившей заявки в соответствующей очереди первой фазы $$(l_{з} +1) $$. Текущее содержимое этой очереди обозначено как $$l_{з}$$.

Блок А8 - запоминание времени постановки заявки в очередь и передача управления блоку А13.

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

Блок А9 - имитация загрузки одного из свободных каналов поступившей заявкой. Очередь свободных каналов уменьшается на $$1 $$ ( $$l_{ск} -1$$ ) . Выбранному каналу присваивается признак занятости - $$1$$. У свободных каналов этот признак равен $$0$$.

Блок А10 - определение и запоминание длительности простоя канала. Накопление времени простоя необходимо, например, для определения коэффициента занятости канала на интервале исследования.

Блок А11 - имитация времени обслуживания заявки.

Случайное время обслуживания $$\tau_{обс}$$ формируется обращением к соответствующему датчику случайных чисел.

Блок А12 - прогноз времени окончания обслуживания $$t + \tau _{обс }$$ и занесение новой пространственно-временной точки в список событий.

Блок А13 - прогноз момента времени поступления очередной заявки. Соответствующий датчик случайных чисел выдает длительность случайного временного интервала между заявками данного типа. Момент поступления очередной заявки вычисляется как $$t + \tau _{инт}$$. Полученная таким образом пространственно-временная точка заносится в список событий.

Управление передается в точку Е - блоку 3 для определения очередного ближайшего события.