Методы статистического моделирования случайных величин по эмпирическим распределениям Текст научной статьи по специальности «Математика»

УДК 621.019

МЕТОДЫ СТАТИСТИЧЕСКОГО МОДЕЛИРОВАНИЯ СЛУЧАЙНЫХ ВЕЛИЧИН ПО ЭМПИРИЧЕСКИМ РАСПРЕДЕЛЕНИЯМ

Э.М. ФАРХАДЗАДЕ, А.З. МУРАДАЛИЕВ, Т.К. РАФИЕВА, У.К. НАСИБОВА АзНИПИИ Энергетики, Азербайджан

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

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

Моделирование непрерывных случайных величин. Обозначим статистическую функцию распределения через Г* (х), где индекс «*» будет характеризовать эмпирический характер распределения случайной величины х. Пусть в результате классификации данных получена выборка из п значений х. Обозначим ее через { х }п. Разместив х в порядке возрастания, придадим наименьшему х первый порядковый номер. Тогда вариационный ряд выборки { х }п будет иметь вид

х 1 < х 2 < х з <.< хп,

где

х 1 = шт{ х }п , хп = шах{ х }п .

Если теперь сопоставить каждому значению х вариационного ряда вероятность Р = 1/п, то аргументу х; статистической функции распределения

будет соответствовать величина Г (х/ ) = ^п. Этот способ определения г1* (х) используется, в частности, при проверке предположения (гипотезы) о соответствии Г* (х) некоторому гипотетическому закону распределения Г (х).

© Э.М. Фархадзаде, А.З. Мурадалиев Т.К. Рафиева, У.К. Насибова Проблемы энергетики, 2008, № 9-10

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

- Г* (х) представляется в кусочно-линейной форме;

- программно моделируется случайная величина \ с равномерным законом распределения в интервале [0,1];

- определяется интервал значений Г* (х), для которого

Г1* (х1 )< \ < Г1* (с+1 ^

- по формуле

(х1+1 - х1 Ж - Г1 (х1)] ( )( Е Л

х = х + г ^—тг = х +(х;+1-х Дп^-1)

[ ([+1)- Г1* (.х1)]

(1)

вычисляется реализация х. Графическая иллюстрация расчета х приведена на рис. 1.

Рис. 1. Графическая иллюстрация расчета случайной величины х

* / \

В работе [1] функция распределения Г (х) задается как

* / V

Г2 (х ) =

0

/ - 1

(х - х1 )

п -1 (п - 1)(х/+1 - х1 )

если х < х1 если х1 < х < хп

если х ^ хп .

(2)

Если преобразовать это выражение относительно х , получим х = х + (х+1 - х )[^(п -1)+(/ +1)].

(з)

Графическая иллюстрация различия этих двух подходов к определению Г * (х) приведена на рис. 2 для реального примера случайных величин длительности простоя первого энергоблока ГРЭС в аварийном ремонте при внезапных отказах за 2002 г. с п = 4 значениями х, равными х1 = 20 ; х2 = 22 ; х 3 = 44 и х 4 = 136.

1

40 80 120 160 ч

* / \ * / \

Рис. 2. Графическая иллюстрация расхождения распределений Г1 (х) и Г2 (х)

Как отмечается в работе [1], существенным недостатком Г2* (х) является то, что моделируемая случайная величина х будет меньше х п и больше х 1. Для

распределения Г (х) реализации случайной величины х удовлетворяют неравенству 0 < х < хп. Но главное, что среднее значение реализаций х не равно

4

среднему значению выборки хСр = ^ х/ /4 . Наши исследования этого вопроса

г = 1

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

С

меньше, чем х ср . Различие тем больше, чем больше неравномерность реализаций х, в частности различие между х п и хп-1, и чем меньше п . Средние значения моделируемых случайных величин х по распределениям Г (х) и Г* (х) соответственно равны х Мр 1 =38 ч. их Мр 2=46 ч. при среднем арифметическом

значении исходных данных х М =56 ч., где индекс «М» обозначает моделирование. Таким образом:

- учет возможности х < х1 приводит к более существенному отличию х с?р и

С

х ср ;

С

- хотя отличие х ср от хср для Г2 (х) меньше, чем для Г1 (х) ,

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

проявления величин х, близких к хп (при хп >> хп-1). Такие расхождения

х п и хп-1 свойственны выборкам, действительное распределение которых асимметрично. В качестве примера асимметрии распределения длительности состояний объектов на рис.3 приведены гистограммы распределения длительности состояний ЭБ 300 МВт.

а)

в)

с)

Рис. 3. Гистограммыраспределения длительностей состояния ЭБ 300 МВт: а) - рабочее состояние; в) - состояние аварийного отключения; с) - состояние резерва

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

*

- при вычислении хп+1 по формуле хп+1 = хср + А • с х,

^ (х хср )

г=1 п -1

возможны два исхода. Если А > 2, что соответствует при экспоненциальном

М

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

х п , и превышение тем больше, чем больше А . Если же А<2, то возможны случаи,

когда расчетное значение хп+1 окажется меньше хп , что приводит к сбою программных расчетов;

- наиболее эффективным оказался подход, при котором хп+1 = х п .

Расхождение между х ^ и х ср оказалось значительно меньше, чем при моделировании по распределениям (х) и Г* (х). В частности, для

рассмотренного выше примера х ^М 3 =58 ч., а расхождение его с хср =56 ч. составило

* / \

не более 5%. Расчет реализации х для Г3 (х) проводится по формуле

х = хг +(хг+! - хг)[(п +1 )£ - /],

где г = 1, (п +1).

Следует отметить, что моделирование реализаций х по Г* (х) и Г2 (х)

М

приводит не только к различию хср и хср , но и к изменению дискретного

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

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

у = 1, тг признаков. Например, признак «узлы трансформатора» включает корпус, вводы, магнитопровод, обмотки, РПН и т.д.

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

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

*

признака П у,/ с I = 1, тг у . Оценку вероятности проявления г -ой РП Qу,г вычислим

*

Несмотря на то, что физическая сущность Гу (V) и известных дискретных

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

- моделируется случайное число § с равномерным распределением в интервале [0,1];

*

- определяется порядковый номер интервала распределения Гу (V), соответствующего вероятности § по формуле

как

а статистическую функцию распределения вероятности

проявления V - ой РП определим по формуле

V

(5)

г=1

где Г* (0) = 0; Г* (тгу )= 1.

Г* (/-1 )< § < Г* (/^

где / = 1, тг,у;

- если условие (6) выполняется, то принимается однозначное соответствие номера г -го интервала и факта проявления г -ой РП.

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

Пусть исходное соотношение вероятностей проявления РП задается рядом

* * А

Q Ц ; Q ! 2 ;..Q • . , каждая из которых вычисляется как Qj; = П /п ; V . Если

J) J * у,^тр у

предположить, что в результате моделирования проявилась V -ая РП, где 1 < V < тг,у, то на очередном шаге моделирования предварительно уточняются

вероятности Q*,г с / =1, тг у по формуле Q*,/ = пу /(п^ -1) для всех значений / за исключением V -го значения:

тг,]

Q = (п у, -1 )(п(,2 -1). Нетрудно заметить, что по прежнему ^ Qj,i = 1.

г=1

Однако «шансы» проявиться у всех признаков, кроме V -го, возросли. На определенном этапе вероятность проявления тех или иных РП становится равной нулю. Ресурс возможности их проявления исчерпан. Процесс моделирования завершается, когда возможность проявления сохранилась лишь у одной РП. Таким образом, вместо случайной стратегии розыгрыша РП вводится стратегия «предрасположения» объекта к наблюдаемому соотношению вероятностей проявления РП. Повышению достоверности этого соотношения способствует алгоритм целесообразности классификацииу-го признака на тгу РП.

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

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

Иллюстрация применения метода «условных вероятностей» приведена в таблице при имитационном моделировании изменения состояний энергоблоков [2]. В заданном интервале времени энергоблок (ЭБ) может находиться в одном из девяти выделенных состояний (девять РП). Два из них (средний и капитальный ремонты) имеют детерминированный (неслучайный) характер.

Вероятность нахождения ЭБ в і -ом состоянии может быть определена как относительная длительность этого состояния и рассчитывается по формуле

пИ

5т а = 2т ї,і,у /АТ,

У=1

где АТ - заданный интервал времени, ч.; т - V -ая реализация длительности

і-ой разновидности у-го признака.

Результаты расчетов 5т у с і =1,9 приведены в первой строке таблицы.

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

5т у,і ^ 5т ],і/(і - 5т у,8 - 5т у,9 ),

где і =1,7; индекс «^» означает соответствие и преследует цель снижения числа переменных и их индексов.

Во второй строке таблицы приведены результаты расчетов значений функции і

распределения Ґ* (5т у )= 2 5т у)Л, с і = 1,7.

V=1

На этом процесс подготовки исходного при моделировании типа состояний

распределения Г* (5т у) завершается. Предположим, что реализация случайного

числа \ оказалась меньше, чем 5туд = 0,72, т.е. первым состоянием ЭБ (с начала

периода моделирования) было рабочее состояние. А результат моделирования длительности этого состояния оказался равным т удд = 500 ч. Относительная

длительность т у 11 равна:

5т у,1,1 = т у,1,1 /[1 - 5ту,8 - 5ту 9 ]8760 = 0,09.

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

5ту, 1 ^ (5туд - 5тудд )/(1 - 5т удд )= (0,72 - 0,09)/0,91 = 0,692, а для остальных состояний - по формуле

5ту,і ^ 5тІї/(1 -5ту,1,1 ).

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

© Проблемы энергетики, 2008, № 9-10

© Проблемы энергетики, 2008, N° 9-10

Таблица

Иллюстрация метода условных вероятностей

Последующее состояние

Предшествующее Рабочее Аварийн. Аварийный ремонт вследствие: Холодный Средний Капит.

состояние, 1 простой отказа при пуске повтори. отказа внезапн. отказа авар. заявки резерв ремонт ремонт I

1 2 3 4 5 6 7 8 9 10

Исходное распределение '/Ы 44,6 0,720 од 0,722 0,8 0,735 0,5 0,743 1,8 0,772 4,9 0,850 9,3 1 - 38 100

Рабочее состояние 5тУ./- % 5тУ./- % '/ы 69,2 0 0,2 0,2 0,007 1,4 0,007 0,9 0,9 0,038 3.2 3.2 0,146 8.7 8.7 0,441 16.5 16.5 1 VI - 500 ч 100 29,5

Состояние холодного резерва '/ы 74.1 74.1 0,873 0,2 0,873 1.5 1.5 0,891 0,9 0,891 3,4 0,891 9.3 9.3 1 10,6 1 ТУ,7-1 = 360 ч 100 84,9

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

5т ji = 5т jil(l - 5т j,1 - 5т j,3).

На основе данных расчета определяется функция распределения F* (5т j ) по формуле

F* -£ 5т Ji .

i—1

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

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

■ проводится коррекция данных 5т j i, приведенных в четвертой строке

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

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

■ определяется функция распределения F (5т j ) (восьмая строка таблицы).

Все последующие расчеты F (5т j ) проводятся аналогично.

Выводы

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

2. Учет особенностей моделирования дискретных состояний объекта достигается разработанным методом «условных вероятностей».

Summary

Statistical modeling of random quantities on empirical distributions concerns to milestones of simulation modeling of process of operation of complex systems. Analysis of existing methods is conducted, their limiting are taped and new methods which implementation confirms boosting adequacy of implementation of flock of simulated random quantities are recommended.

Литература

1. Аверилл М.Лоу, Кельтон В.Дэвид Имитационное моделирование. Классика CS / Перевод с англ. 3-е изд. - Питер; Киев: Издат. Группа, 2005. - 847 с.

2. Фархадзаде Э.М., Мурадалиев А.З., Рафиева Т.К. Имитационное моделирование состояний энергоблоков. - Баку: ЭЛМ, 2005. - С. 23-32

Поступила 28.02.2008