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

УДК 519.87

Муравьев1 И.И., Острейковский2 В.А., Шевченко2 Е.Н.

1ПУ «СургутАСУнефть» ОАО «Сургутнефтегаз», Сургут, Россия 2«Сургутский государственный университет», Сургут, Россия

МОДЕЛИ ОЦЕНКИ ФАКТОРА ВРЕМЕНИ В ТЕОРИИ ТЕХНОГЕННОГО РИСКА ДИНАМИЧЕСКИХ СИСТЕМ

*Работа поддержана РФФИ (проект 14-01-00230)

Введение

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

Постановка задачи

В работах [1-3,7] разработана классификация моделей техногенного риска сложных динамических систем, в которых детально сформулирован теоретико-множественный подход к анализу характеристик множеств , и при математическом моделировании техногенного риска. В общем виде графическая интерпретация соотношения множеств риска , вероятностей исходных событий (отказов, аварий, катастроф) 2, ущерба от них времени Т приведено на рис.1

во

Рисунок 1 - Множества К, С, 2 и Т в моделях техногенного риска при количественной оценке безопасности сложных систем

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

к = Н{2*С*Т}, (1)

где: реализуют следующие отобра-

жения: Н : 2 * С-> К ± при независимых между собой 2 и С; Н2:2 * С/(2 К2 при коррелированных 2 и С; Я3:2*С*Г->К

Рисунок 2 - Ущерб как случайный процесс

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

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

дующими свойствами:

процесс изменения во времени рас-

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

длительность наблюдения за системой в течение конечна;

энергетический спектр случайного процесса сплошной и отличен от нуля на всей оси частот ;

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

?? = / (Т) ограничен, причем т0 макс«Т. Под интервалом корреляции принимается промежуток времени, за который полностью затухают корреляционные связи между отдельными составляющими К ( {), 2 (0 и С ({). При нестационарном характере изменения случайной функции интервал корреляции в общем случае зависит от рассматриваемого момента времени, т.е. т0 = / ({) .

Нестационарные случайные процессы изменения целесообразно представлять в виде суммы нескольких процессов К = / (Т) :

, (2) где А ({) - нестационарный случайный процесс, характеризующий необратимые изменения в системе в результате старения, изнашивания, регулирования; В ({) - стационарный случайный процесс, характеризующий обратимые изменения из-за колебаний внешних условий при эксплуатации системы; - стационарный случайный процесс ошибок измерений 2 ( 0 и С ( {).

Рассмотрим методы оценки фактора времени, приведенные в [4], при планировании многофакторных испытаний на надежность сложных критичных систем (параметра 2)-Основная часть

Метод фактора времени при моделировании риска в условиях ортогонального дрейфа.

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

? ^ = М[%] = У ( 0 (2)

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

по выбранной системе ортогональных функций, часть их используется для описания дрейфа, а часть — для варьирования управляемыми факторами, т.е. образования плана. Предложенный Боксом [4] метод построения планов для оценки поверхности отклика в условиях дрейфа основан на использовании полиномов Чебышева, Сущность метода заключается в следующем. Функция дрейфа ? = У({) на интервале описывается полиномом порядка :

(3)

где — ортогональные полиномы Чебышева, определенные для равностоящих абсцисс наблюдений интервала Т.

Чтобы выделить влияние управляемых переменных ( и С независимо от дрейфа, вектор-столбцы планирования формируют в виде линейных комбинаций всех оставшихся полиномов порядка к > d. В этом случае модель имеет вид

Ш С, Т) = И(Х, Т) = а0Ь0(Т) + а1Ь1(Т) + - + ас1Ьс1(Т) +

г>о(0

Ь1х1 + Ь2х2-\------Ькхк, (4)

где X С (хС; Ъ0(¡) определяет смещение поверхности отклика во времени, а коэффициенты Ь1 — ее наклон в исследуемой области К ( (},С). Анализ уравнения (4) сводится к оценке значимости эф-

фектов дрейфа и управляемых факторов по ¥-критерию. Каталог планов для наиболее употребительных значений Ыг dr к приведен в [4].

Применение ортогональных полиномов Чебышева особенно эффективно, когда отсчеты времени делаются через равные интервалы ДС, а для каждого значения t ставится одно и то же число параллельных опытов. Шкала времени при этом преобразуется в последовательность натуральных чисел 1, N. В этом случае ортогональные полиномы равны [4]:

<р0 = 1;

<Р2 =№[(£-О2-(«2- 1)/12]; <р3 = - О3 - (ЗМ2 - 7)({ - 0/20]; <г>4 = - О4 - (З«2 - 13)({ - 02/14 + 3(«2 - 1)(«2 - 9)/560]; <р5 = - О5 - 5(7«2 - 7)({ - 02/18 + (15Л[4 - 230«2 + 407)0 - 0/1008];

(5)

Коэффициенты зависят от числа наблюдений N. Полиномы (5) табулированы [4].

Недостатки метода Бокса:

1) планы, основанные на применении ортогональных полиномов Чебышева, могут быть использованы только в том случае, когда отсутствуют корреляционные связи типа «время—фактор»;

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

3) необходимо, чтобы на стадии предварительных исследований был выявлен порядок полинома дрейфа dr что само по себе эквивалентно проведению громоздкого эксперимента.

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

Функцию отклика К (((,С,Т) системы в условиях динамической задачи можно представить в виде:

Шс.Т) = яМХЛ + ъиьМх^^ь^х^ + 1 ь1, ( о х2 + - ■ ■■, (6)

ХС (¿*С.

Целесообразно преобразовать выражение (6) в уравнение с числовыми коэффициентами регрессии. Принципиально это можно сделать несколькими путями.

Рассмотрим один из методов для случая линейной модели:

Й ((2, с, т) = к 0 ((2, с, Т) + 1Ь1 (0 х . (7)

Представим К 0 ( (,С,Т) и Ь1 ( 0 в виде функциональных рядов

(8)

= ^с0ц<рц(0,

_ и=1

где <ри&), и = 1, N - коэффициенты функции; с0и,с0и,и = Х N,1 = %1 - постоянные коэффициенты. Тогда

К ( (,С,Т)=£ ^1Сои<рМ + 1к=1'Еи=1С1и<Ри(Ох. (9)

Проинтегрируем выражение (9) справа и слева в пределах (0;7):

I Й(<2,С,Г]№ = I I <Рц(0^Х;.

^ и = 1 ^ 1 = 1 и = 1 ^

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

11((2,с,т)(И = с*0и = с0и ск =

Jo Jo Jo

— С0и — С0-Jо

Тогда получим выражение

-II = 1. £-41 = 1

(10)

В этом уравнении коэффициентами регрессии являются величины и ,

Оценки этих коэффициентов определяются соотношениями:

(11)

Функции и значение выбираются такими,

чтобы

¡0 <Ри(0& Ф 0, и =

Модели оценки фактора времени в эксперименте, когда одна из контролируемых переменных -время

Представим выражение для функции отклика в виде [4]:

К (((, С, Т) = [в (((, С, в) + V (((, С) ] * Г (с) + £ ( 0 , (12)

где и -

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

мк<?,с)] = о,м[£(0] = о л

где и - между собой независимы;

КЕ(^{2) - корреляционная функция риска.

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

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

Случай 1.

Кк (^Л2) «г(С)(((,С) Г'(1)ЛЕ[0,Т (14)

или

(15)

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

Имея реализацию К [ ((,С) ) можно определить значения при заданных по

формуле:

e[(Q,c)t,в] + vt = [/07(огак]-1 /отк[«?, о„тт,

(16)

минимизируя интеграл

[ {№(<?, с)ов]+ ц]'№- я[(<г,с)о 0}2^ ь

При конечном числе замеров реализации вместо выражения (16) имеем:

9№С)1,в] + и( =

(17)

Количество замеров и положение точек должны быть таковы, чтобы матрица

а-

/о <РиШ(

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

и задача становится эквивалентной следующей регрессионной задаче:

Матрица 0„ ( 2, С) предполагается известной. На основании измерений случайной величины получаются оценки поверхностей или оценки параметров В.

Если зависимость от неизвестных параметров линейна то для оценок пара-

метров регрессии можно построить наилучшие линейные оценки:

В = М" 1 IV, (18)

где

к 1 = 1

к

ш = ^ртс)1]й-1тсж-, Л(<2, С)(] = \\<рШ С)г].....<ртМ. С)Л1-

Случай 2. Неравенство (14) не выполняется, но корреляционная функция ( {2) известна. В этом случае моделирование для регрессионной задачи (12) так же удается свести к моделированию в случае одновременного измерения нескольких величин [4].

При известной функции наилучшие ли-

нейные оценки для величин ^ = 0 [ (2, С) ¡, В 1 + 1 равны:

z, = iir'Wi,

где

(19)

(20)

wi = ZTj с)<- t/<];J

Причем дисперсионная матрица оценок Z¡ при заданном ^¡обратна матрице А':

б (Z¡/ v¡)=A Г 1 ■ (21)

Из выражения (20) следует, что

(22)

При этом D (fí¡) = А'-1, y¡ и fí¡ независимы, и M[Z] = 0[(<?,С%В]; D[Zt] = Z)„[((?, C)¡] + Af1.

Таким образом, задача сведена к случаю 1 лишь после произведения замены D„ [(<2, С)] на

D„ [( 2 ,С) ] + АТ Ч

Случай 3. Неравенство (14) не выполняется, корреляционная функция ( t^ t2) неизвестна. В этом случае непосредственное применение формулы (19) невозможно. Однако во многих случаях имеется информация о скорости убывания | (t^ t2) |. Если измерения проводятся через промежутки времени, за которые становится близким к нулю, а амплитуда внешних помех не зависит от t, т.е. (t^ t2) = 0"| (это требование обычно выполняется при не слишком больших интервалах [0,Г] ) , то выражение для А 1 и Напринимают вид:

А1*<т-2Xf= i/ (í)/ '(í) ; (23)

Hi * <7-25У=1/(^)Я[(2,C)¡,t] ; (24)

Если значение неизвестно, то в выражении (23), (24) следует подставлять его оценку

s2 _ Ijli[«[(O.C)¡.tj]-Z'(g,C)¡/(tj)] m-n

(25)

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

, можно воспользоваться оценкой для :

Г1

(г, -1)0(2,) = - гп)(2п - 2Г1)':

1=1

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

данном ( 2, С) ¡.

Точку (2, С) ¡стремятся выбрать так, чтобы она принадлежала оптимальному плану.

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

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

лишь знание дисперсионной матрицы D(2t/ v¡), соответствующей использованному методу [4].

В дополнение к изложенному в п.1-3 могут быть использованы и другие методы количественной оценки фактора времени. В частности, в [57] были эффективно применены к поставленной в данной статье задаче методы спектрально-корреляционной теории случайных процессов регрессионного анализа и теории катастроф.

Метод оценки фактора времени с применением спектрально-корреляционной теории случайных процессов

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

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

Рассмотрим случай с линейной дискретной моделью риска R (t), который задается Фурье-разложением по конечному числу частот:

к (t) = YljL i ( О/ С О S t + bbj si n t) , (26)

где и - независимые центрированные, нормально распределенные величины с дисперсией af = M[aj] = M[bj\ которая является функцией частоты ш. Корреляционная функция случайного процесса равна

М [R (t)R (t + т)] = tfR (т) = Sf=1 с/2 cos ш/т. (27) Согласно теоремы Хинчина-Винера [8] спектр случайного процесса и корреляционная функция являются сопряженными Фурье-функциями. Поэтому

Sf = ± [J?R(0) + 2 ЯУ KR(t)KR(t + г)},

и Zfp

(28)

выборочные оценки. Для выборки

объема N

КлСО = —^ÍKrWRÜ + т),

(29)

где т принимает лишь дискретные значения.

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

KR (0) + 2 ^ ЛТКЕ (т) eos (üjT

где

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

положений об исследуемом случайном процессе.

Для нормального случайного стационарного процесса средняя квадратичная погрешность вычислений корреляционной функции, возникающая из-за конечной длины реализации, при условии М[К(С)1 = = 0 равна [9]

< =^/вГ~Т(Т ~ т ~ е)№(е) + Кн(е + т) + кв(в -т) ] й ( 0 ) , (30)

где

где КЕ — корреляционная функция ВПО К ( С); Т — длина реализации (время наблюдения); — временной сдвиг; в — приращение временного сдвига.

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

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

(32)

в [4] для значений получена следую-

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

о1„<гА—„, (33)

приближенная зависимость (34) имеет вид:

икя - 'яг'

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

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

= D В

■ ß, KR (т) = D (а е

-с|т| _ hp~d 1Т1

). (34)

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

\Кц(ттах) < 0,051/Сд (0) 11.

При колебательном характере «хвоста» кривой корреляционной функции вычисления прекращаются в точке отстоящей по оси абсцисс на рас-

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

вида

КЕ (т)= В(^е-сМ-1е-лМ ) (35)

где - параметр формулы аппроксимации, характеризующий наклон функции.

Введя понятие относительной средней квадратичной погрешности определения дисперсии г]в = ав/0 и подставив его в формулу (36), получим

Т = — . (37)

СЧО

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

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

Т = —

(38)

V

у = тс; Рв = —е~у = (1 — е~г)2 .Требуемое

где

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

Заключение

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

Установить равную для всех режимов функционирования системы продолжительность наблюдения;

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

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

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

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

ЛИТЕРАТУРА

1. Королев В.Ю., Бенинг В.Е., Шоргин С.Я. Математические основы теории риска: Учеб. Пособ. -М.: ФИЗМАТЛИТ, 2007. - 544 с.

2. Острейковский В.А. Теория техногенного риска: математические методы и модели: монография, Сургут, гос. ун-т ХМАО - Югры. - Сургут : ИЦ СурГУ, 2013. - 320 с.

3. Острейковский В.А., Шевченко Е.Н., Микшина В.С. Количественная оценка риска в техногенной безопасности сложных динамический систем: монография // Итоги науки. Том 1. - Избранные труды международного симпозиума по функциональным и прикладным проблемам науки. М.: РАН, 2013. С. 12 -31.

4. Острейковский В.А. Многофакторные испытания на надежность. - М.: Энегрия, 1978. - 152 с.

5. Острейковский В.А. Теория надежности: Учебник для вузов - 2-е изд., испр. - М. : Высш. шк., 2008. - 463 с.: ил.

6. Жакот А.Д., Острейковский В.А., Челнокова Е.В. Математические модели регрессионного анализа и теории катастроф синдрома дыхательных путей / Надежность и качество. Тр. междунар. симпозиума / под ред. Н.К.Юркова. - Пенза : Изд-во ПГУ, 2005. - с.467 - 472.

7. Острейковский В.А. О некоторых классах моделей риска в теории техногенной безопасности // Надежность и качество. Тр. междунар. симпозиума: в 2 т. / под ред. Н.К.Юркова. - Пенза : Изд-во ПГУ, 2013. - 1 т. - с.46 - 49.

8. Гихман И.И., Скороход А.В. Введение в теорию случайных процессов. Изд. 2-е, Главная ред. физ.-мат. лит. изд-ва «Наука», М., 1977. - 568 с.

9. Пугачев, В.С. Теория случайных функций и ее применение к задачам автоматического уравнения. - Изд. 3-е, испр. - М. : Гос. изд-во физ.-мат. лит., 1962. - 883 с.

10. Артемов И.И. Исследование влияния дефектной структуры материала болтового соединения на процесс ослабления затяжки / Артемов И.И., Кревчик В.Д., Суменков С.В. // Новые промышленные технологии. 2002. № 5-6. С. 67.