5. ПРИМЕРЫ МОДЕЛИРОВАНИЯ СЛУЧАЙНЫХ ПРОЦЕССОВ МЕТОДОМ МОНТЕ-КАРЛО

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

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

Требуется найти следующие характеристики работы устройства:

— вероятность исправной работы устройства как функцию времени;

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

— среднее время работы системы, т. е. среднее время, которое устройство будет проводить в работающем состоянии;

— среднее число запасных узлов которое будет израсходовано, а также среднее число у, израсходованных запасных узлов

Рис. 8.16

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

берем случайное число R от 0 до 1 и подвергаем его преобразованию (рис. 8.16)

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

Рассмотрим четыре параллельные оси с одним отсчетом времени (рис. 8 17) На оси (1) мы будем отмечать состояние первого узла (жирной линией — «работает», тонкой — «отказал»). На оси (2) также отмечаются состояния второго узла, на оси (3) — третьего, на оси (4) — состояние системы в целом («работает», «не работает»).

Так как момент отказа узла нам известен, то мы можем сразу же заполнить ось (3). После этого будем заполнять (1) и (2). Сначала разыграем время в течение которого будет работать основной узел — для этого мы воспользуемся функцией распределения

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

Если в момент окончания этой замены третий узел еще работает то снова разыгрываем значение времени работы первого запасного узла Г, с помощью функции (5.1) и, после этого — опять время замены этого узла с помощью функции (5.2).

Предположим, что момент окончания этой замены (как показано на рис 8.17) оказался после момента

Тогда при розыгрыше нового значения времени безотказной работы узла мы должны будем воспользоваться уже не функцией (5.1) а новой условной функцией распределения

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

Рис. 8.17

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

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

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

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

Полученное таким образом значение нужно прибавить к уже прошедшему времени

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

Рис. 8.18

Рис. 8.19

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

Таким образом, разыграна одна реализация нашего случайного процесса. Разумеется, если моделирование производится на ЭЦВМ, никаких графиков, осей и участков строить не нужно; розыгрыш обеспечивается приведением в действие машинного расчетного алгоритма, который сочетает единичные жребии — розыгрыши со сравнением между собой моментов осуществления разных событий (что произошло раньше — восстановление первого (второго) узла или выход из строя третьего?).

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

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

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

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

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

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

где — число запасных узлов израсходованных в реализации.

Аналогично определяется

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

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

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

— среднюю долю пораженной площади цели:

— дисперсию доли пораженной площади цели:

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

— математическое ожидание числа ракет, которые не причинили цели никакого ущерба (попали мимо).

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

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

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

Предположим, что моделирование производится на ЭЦВМ. Тогда наиболее удобным способом розыгрыша пары нормальных величин будет описанное в § 3 сложение нескольких независимых случайных чисел от 0 до 1 с последующей перенормировкой. При таком способе координаты точки попадания могут быть разыграны по формулам:

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

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

Рис. 8.20

Рис. 8.21

Деля полученное в реализации значение на площадь цели, получим долю пораженной площади в данной реализации:

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

Средняя доля пораженной площади:

Дисперсия доли пораженной площади:

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

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

где — число ракет, не причинивших ущерба цели в реализации.