Модели и методы. Модели и методы построения вероятностно статистических оценок для мониторинга показателей надёжности в диспетчерском управлении транспортом газа
 Скачать 3.85 Mb.
|
|
8,54 12,94 33,86 42,98 37,10 1,6 4,06 6,39 10,32 27,63 34,18 31,46 1,7 2,90 4,85 8,32 22,45 27,53 26,91 1,8 2,11 3,74 6,78 18,22 22,42 23,20 1,9 1,56 2,91 5,58 14,80 18,45 20,14 2 1,17 2,30 4,63 12,06 15,32 17,60 2,1 0,88 1,83 3,87 9,86 12,83 15,47 2,2 0,68 1,47 3,26 8,10 10,82 13,67 2,5 0,33 0,81 2,03 4,63 6,76 9,70 2,6 0,26 0,67 1,75 3,88 5,85 8,73 2,7 0,21 0,56 1,52 3,27 5,09 7,88 2,8 0,17 0,47 1,32 2,77 4,44 7,14 2,9 0,14 0,40 1,16 2,36 3,90 6,49 3 0,11 0,34 1,02 2,02 3,44 5,91 3,1 0,09 0,29 0,90 1,73 3,04 5,41 3,2 0,08 0,25 0,80 1,50 2,70 4,96 Относительная погрешность формулы (19) относительно параметра не превышает 10% при 1,5 ; формулы (21) для V C – не превышает 10% при 2, что вполне допустимо при инженерных расчетах. 60 Абсолютная погрешность формулы (24) для относительно параметра не превышает 3% при 1 иона уже меньше 1,5% при 2 (см. Рисунок 8). Рисунок 8. Графическая зависимость абсолютной погрешности формулы (24) Приближенные формулы ускоренных расчётов дисперсии и коэффициента вариации, приемлемые для экспресс-анализа С помощью компьютерного моделирования в профессиональном математическом пакете Mathematica были получены следующие приближенные формулы для ускоренных расчетов дисперсии и коэффициента вариации 2 2,345 1 1 D , 0,934 С (25) Абсолютная погрешность приближения (25) для D относительно параметра меньше 3% при 1,5 ; формулы для Сне превышает 0,7% при 1 (см. Рисунок 9). Данные соотношения могут пригодиться при инженерном экспресс– анализе. Рисунок 9. Графическая зависимость абсолютной погрешности Сот 2 3 4 5 6 7 0.001 0.002 0.003 0.004 0.005 0.006 0.007 61 Кроме того, в результате графического моделирования в пакете Mathematica позволило получить для ( , ) GW следующую приближенную формулу 2,86 ( 1) 2 S V A C − + , (27) Нов силу (25) имеем 0.934 2,86 0,86 S A − − (28) Рассмотрим некоторые частные случаи. В случае показательного распределения 1 = , 1 V C = и 2 S A = , что вполне согласуется. Распределение Рэлея ( 2 = ), тогда из теории [46]: 4 0,523 V C − = и 3 2 ( 3) 0, 631 (4 ) S A − = − , а по формуле (28) 0,523 V C и 0, Если же 3, 6 = , когда 2 ( , ) ( ; ) GW N a (о чем речь впереди, см. п. 2.5), то 0 0,30069 S V A C = , что тоже вполне согласуется с классической теорией случайных величин (см. Таблицу 2). Таблица 2. Диапазон коэффициентов вариации для стандартных распределений Пределы изменения коэффициента вариации V C Закон распределения СВ 0 Нормальный 0,3 0, 4 V C Гамма–распределение 0, 4 1 V C Гнеденко-Вейбулла Экспоненциальный Указанную зависимость (27) представлена графически на Рисунке 10. 62 Рисунок 10. Графическая зависимость асимметрии от коэффициента вариации Показатели надёжности невосстанавливаемых объектов в модели Гнеденко-Вейбулла распределения отказов Исходным показателем надёжности невосстанавливаемых технологических объектов ГТС, оценка которого производится на основе обработки собираемых статистических данных, является интенсивность отказов ( ) t , определяемая следующим образом ( ) ( ) 1 ( ) f t t F Отметим преимущества введенной величины − по известной интенсивности ( ) t несложно оценить остальные показатели надёжности, те. функция интенсивности отказов является одной из форм закона распределения случайной величины 0 0 0 ( ) ( ) ( ) exp ( ) , ( ) 1 ( ) exp ( ) t t t f t P t t s ds P t f s ds s ds = − = − = − = − – функция ( ) t наглядно описывает все этапы жизненного цикла функционирования объекта (см. рис, в частности кривая интенсивности отказов достаточно просто характеризует процессы деградации 5 10 15 20 25 30 20 40 60 80 63 − интенсивность отказов может быть легко определена экспериментальным путем. Для распределения Гнеденко-Вейбулла функция интенсивности отказов имеет вид 1 ( ) t t − = (29) Распределение Гнеденко-Вейбулла позволяет аппроксимировать экспериментальную кривую интенсивности отказов на каждом из основных периодов функционирования системы [3]. В частности, период приработки отвечает распределению Гнеденко-Вейбулла с параметром (0;1) ; период нормальной эксплуатации – с параметром 1 и период старения – с параметром 2 . Это следует из исследования функциональной зависимости от параметров и интенсивности отказов ( ) t (см. Рисунок 11). Рисунок 11. Аппроксимация кривой интенсивности отказов распределением Гнеденко-Вейбулла (пунктирная линия) Применимость распределения Гнеденко-Вейбулла в описании функционирования объектов с позиций жизненного цикла, в том числе процессов деградации Здесь следует еще раз обратить внимание на тот факт, что в модели отказов по распределению Гнеденко-Вейбулла, для участка устойчивого функционирования, коэффициента для третьего участка коэффициент 64 больше 2 (см. Рисунок 11). Разрыв между численными значениями коэффициента, очевидно, объясняется тем, что, практически, деградация как эволюционный процесс является непрерывным, а его возникновение происходит не дискретно, а непрерывно. Таким образом уже на участке нормальной работы ближе к началу третьего этапа, начинается процесс деградации. Из методики формирования оценок надёжности оборудования следует, что данное значения параметра характеризует нахождения оборудования в конце II- этапа эксплуатации, на котором этап нормальной эксплуатации переходит в преддеградационное состояние (Рисунок 12). Рисунок 12. Интенсивность отказов оборудования и граничные значения интервалов для параметра . Главной проблемой в решении подобных задач является определение критического значения crit t – при котором можно считать, что оборудование перешло из периода нормально эксплуатации в период деградации или старения. Из теории известно, что при 2 наступает деградационный период, нона практике такого четкого разграничения на принадлежность реального состояния оборудования к одному из трех периодов эксплуатации не бывает. Характер 65 перехода от первого этапа ко второму, также, как и от второго к третьему - достаточно плавный, например, при значениях близких к 2, сложно идентифицировать состояние оборудования по оценки показателей надёжности. Можно наблюдать, что в конце второго периода уже начинают сокращаться интервалы между отказами, что говорит о начале третьего периода. Поэтому, необходимо применять комплексный подход к оценке показателей надёжности с использованием других методов, например, нечеткой логики. Прогнозирование наступления момента очередного отказа элементов системы ГПА Моменты начала деградационных процессов при эксплуатировании объектов ГТС, в частности ГПА и САУ ГПА, можно найти на основании изучения теоретической плотности вероятности отказа оборудования, подчиненной двухпараметрическому распределению Гнеденко-Вейбулла ( , ) GW Сделаем эвристическое допущение о том, что максимальный рост вероятности отказа оборудования, который соответствует первой точке перегиба функции плотности вероятности отказа, находящейся левее математического ожидания случайной величины , соответствует моменту времени следующего отказа. Проведя несложные преобразования, получим формулу для нахождения точки максимального изменения плотности вероятности отказа, соответствующей первой точке перегиба плотности распределения (2) ( , ) GW : 1 3( 1) ( 1)(5 1) 2 crit t − − − − = . (30) Данное выражение определено при значениях 2 . При единичном значении параметра масштаба ( 1 = ) график crit t как функции от параметра формы имеет вид 66 Рисунок 13. График зависимости crit t от параметра при 1 = . Отметим, что в нашем случае t M crit . Несложно убедиться в том, что 1 crit t − , (31) ввиду того, что 1 3( 1) ( 1)(5 1) lim 1 2 →+ − Следовательно, в терминах интенсивности отказов, получаем 1 2 ( ) crit t − (32) Итак, пусть имеются эксплуатационные данные об отказах активных элементов ГТС , 1,.., k k n = . Формируем последовательность , 1,.., 1 k t k n = − где 1 k k t k + = − . Предполагая, что распределение времени между отказами подчиняется двухпараметрическому распределению Гнеденко-Вейбулла с функцией распределения (1) найдем точечные оценки неизвестных параметров и . Среднее время безотказной работы, те. наработка до отказа, составляет 1 1 1 1 T M = = + . Следовательно, значения должны k t лежать в окрестности 1 T в период нормальной эксплуатации, характеризующийся постоянством 4 6 8 10 0.9 1.0 1.1 1.2 1.3 1.4 67 интенсивности отказов. Процесс старения (деградации) означает нарастание интенсивности отказов, интервалы между отказами уменьшаются, что соответствует приближению k t к Детально изучим первую точку перегиба распределения Гнеденко-Вейбулла с помощью асимптотических методов. Справедливо следующее разложение при больших значениях параметра формы , полагая параметр масштаба 1 = : 1 1 3( 1) ( 1)(5 1) 1 3 1 1 3 1 5 2 2 crit t − − − − = = − − − − = ( ) 2 2 2 1 1 3 ln 3 5 2 3 5 1 1 ln 2 3 5 1 6 1 ln ln 2 ln 3 5 2 2 5 3 e − − − + + − + − + + − − ( ) 2 0,962 2, 688 1 − + → (33) в силу того, что 1 1 1 ( ) 0 1 5 5 → → − − → 2.3. Асимптотическое исследование функции средней остаточной наработки, дисперсии и коэффициента вариации остаточной наработки Средняя остаточная наработка (остаточное время жизни) как функция от времени является важной характеристикой (или даже мерой) процессов старения в приложениях теории надёжности. Ее используют в актуарной математике и страховании, в медицине и биологии. Теоретические свойства ее были 68 рассмотрены Коксом в 1962 [60]. Обзор по теории и приложениям средней остаточной наработки имеется в книгах [90, 133]. Остаточная наработка и функция средней остаточной наработки Пусть время T безотказной работы элемента является случайной величиной, подчиненной двухпараметрическому закону распределения Гнеденко-Вейбулла ( , ) GW с функцией распределения (1), плотностью (2) и математическим ожиданием ( ) 0 0 1 Рассмотрим условную случайную величину ( ) t X T t T t = − , которая называется остаточной наработкой или остаточное время безотказной работы, являющейся интересным объектом для анализа как в теории надёжности, таки в проблемах качества и безопасности. Найдем математическое ожидание указанной случайной величины, те. функцию средней остаточной наработки или математическое ожидание остатка долговечности) на отказ. Имеем 1 ( ) ( ) ( ) ( ) t X F t P t T t F P T t T t − + = − = = 1 ( ) 1 ( ) ( ) ( ) ( ) 1 ( ) ( ) 1 1 , ( ) F t F t P t F t F t F t F t P t − − + + − − + + − = = = где ( ) 1 ( ) P x F x = − обозначает вероятность безотказной работы. Тогда 0 0 0 0 1 1 ( ) ( ) ( ) ( ) ( ) ( ) t t X MX x dF x x dP t x x P t x P t x dx P t P t + + + + = = − + = − + + + = 0 1 1 ( ) ( ) ( ) , ( ) ( ) t P t x d x t P y dy P t P t + + = + + так как, по правилу Бернулли-Лопиталя: 2 2 1 ( ) ( ) lim ( ) lim lim lim ( ) 0, 1 1 x x x x F t x f t x x P t x x f t x x x →+ →+ →+ →+ − + − + + = = = + = − 69 в предположении существования второго момента распределения, что, очевидно, верно, для распределения Гнеденко-Вейбулла. Итак, функция средней остаточной наработки определяется следующим образом ( ) 1 1 ( ) ( ) ( ) 1 ( ) ( ) 1 ( ) t M T t T t P x dx F x dx P t F t t t + + = − = = − − . (34) В случае существования плотности распределения можно получить другое выражение для нахождения ( ) t : ( ) ( ) ( ) ( ) x f x dx t t M T t T t t P t + = − Выразим через ( ) t базовые показатели надёжности, а именно интенсивность отказов ( ) t и вероятность безотказной работы ( ) P t , и как следствие, функцию плотности времени безотказной работы ( ) f t . Из соотношения (34) следует, что ( ) ( ) ( ) 1 ( ) 1 ( ) t F t F x После дифференцирования получаем ( ) ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) 1 ( ) 1 ( ) t t F t t f t F x F t + − + − = − = − Следовательно, 1 ( ) ( ) ( ) t t t + = , и, тогда, по известному в теории надёжности соотношению [78]: 0 ( ) ( ) t d P t e вероятность безотказной работы может быть выражена через функцию ( следующим образом 70 0 0 0 0 0 0 ( ) 1 1 ( ) 1 ln ( ) ( ) ( ) ( ) ( ) ( ) (0) ( ) ( ) t t t t t t d d d d d t P t e e e e + − + − + − − = = = = , то есть, 0 ( ) (0) ( ) 1 ( ) ( ) t dz z P t F t e t − = Дифференцированием последнего соотношения получаем, что функция плотности времени безотказной работы элемента имеет вид 1 0 ( ) 2 (0) ( ( ) 1) ( ) ( ) t z dz t f Применяя правило Бернулли–Лопиталя к выражению для ( ) t , получаем соотношение ( ) ( ) 1 ( ) ( ) 1 ( ) 1 ( ) 1 1 lim ( ) lim lim lim lim ( ) ( ) 1 ( ) t t t F t f t F x dx F t t t f t t t t F t →+ →+ →+ − − + − − − = = = = → + → + − , устанавливающее интересную связь между предельными значениями ( ) t и ( Отсюда следует предельное соотношение для средней остаточной наработки в случае распределения Гнеденко-Вейбулла ( , ) GW : , 0 1 1 lim ( ) , 1 . 0, 1 t t →+ + = = (35) Отметим, что несмещенной статистической оценкой ( ) t для выборки 1 2 , ,..., n t t t является эмпирическая функция [90]: 71 ( ) ( ) 1 1 ( ) [ ] ˆ ( ) 1 1 ( ) [ ] n i i i n n n i i t t t t t P t t t = = − = − где [ ] i t t – индикаторная функция 1, [ ] 0, i i i t t t t t t = Асимптотическое разложение средней остаточной наработки Теперь найдем аналитическое представление для ( ). t Имеем, исходя из определения ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 ( ) 1 ( ) ( ) 1 1 ( ) 1 0 0 0 t t t F t F x dx x t t e e dx t t t x x x e e dx e dx e e dx − + − + − = Для вычисления интеграла используем неполную гамма-функцию (см. 8.350(1) из [33]): ( ) 1 , 0 x t a a x e t dt − − = (36) По формуле 3.381(8) из [33] имеем 1 1 ( ) , (Отсюда 72 ( ) 1 1 1 1 ( ) ( ) 1 , t t e t = + − (37) Преобразуем (37), используя свойства гамма-функции Эйлера ( ) ( ) ( ) ( ) ( ) ( ) 1 1 1 1 ( ) 1 , 1 1 1 1 1 1 1 1 , , t t t t e t e t e t = + − = = − = Отсюда ( ) ( ) 1 1 ( ) , t t e t = , (38) где ( ) 1 , a y z a z y e dy + − − = – неполная гамма-функция (см. 8.350(2) из [33]). Получим еще одно аналитическое представление средней остаточной наработки через функцию Куммера-Похгаммера, воспользовавшись аппаратом специальных функций. Справедливо следующее соотношение 1 ( ) ( ) 1; 1; ( ) 1 1 0 t x t e dx te F t − − = + , где ( ) 1 ( ) ( ) ; ; 1 1 1 (1 )( ) ( ) ( ) ! 1 0 0 k l k k k F x x x l l k k k k l − + + + + = + = + + + стандартное обозначение вырожденной гипергеометрической функции города или функции Куммера-Похгаммера, относящейся к классу специальных функций [6, 74]. Доказательство приведено в Приложении Б. Следовательно, 73 1 1 1 ( ) ( ) ( ) 1 1; 1; ( ) 1 1 t t t e t e F t − = + − + (40) Получим выражение для ( ) t в виде обобщенного степенного ряда. Применяя формулу (39) получаем 1 1 ( ) ( ) 1, 1; ( ) 1 1 1 1 0 1 k t t F t e k k + − + = + = + +Таким образом, 1 1 ( ) ( ) ( ) ( ) 1 1 1 0 0 0 1 1 k k k t t t x t t e dx e te k k k k + + + − − − = + = = = + +где ( ) k a обозначает символ Похгаммера [6]: ( ) ( ) ( ) a k a k a +Окончательно приходим к следующему разложению ( ) ! ( ) 1 0 ! 1 0 1 k t t k t T k k k + = − = + + , (41) где 0 T определено ранее. Приданное соотношение дает 1 ( ) t const = = , соответствующее показательному распределению. Рассмотрим важный частый случай распределения Гнеденко-Вейбулла, когда значение параметра формы 2 = , так называемое распределение Рэлея. 74 1 2 2 1 2 ( ) 1 2 0 t t x t e e dx − − = + − Так как 1 1 2 2 + = , ( ) ( ) ( ) 2 2 2 1 2 2 2 0 2 0 0 x t t x e dx e d x t − − = = , где ( ) 2 2 1 0 2 0 x e dx − = – функция Лапласа. Следовательно, ( ) ( ) 2 2 1 1 ( ) 2 2 0 0 2 2 t t t e t e t = − = − Данную зависимость представим в графическом виде (см. Рисунок 14), считая, без ограничения общности, что значение параметра масштаба 1 = . Рисунок 14. Функция средней остаточной наработки для распределения Рэлея. Укажем оценку погрешности в случае приближенного вычисления ( ) t по найденному аналитическому разложению. Имеет место следующее утверждение 1 2 3 4 5 0.2 0.4 0.6 0.8 75 1 0 1 1 0 ( ) ! ( ) 1 ( ) ! k N N k k t t k t T R t k − + + = = − + , где 1 1 ( ) 1 1 ( ) 1 ! 1 ( ) N N t N R t t N N t − + + + + − , 1 ( ) N t + Замечание Очевидно, что 1 lim ( ) 0 N N R t − →+ = при фиксированных , , t Доказательство приведено в Приложении Б. Изучим асимптотическое поведение ( ) t при t → + . Справедливо асимптотическое представление для ( ) t : 1 1 1 ( ) 1 ( 1) 1 ( ) , ( ) 1 k k t t t t k k + − − + − − → + = (42) Или более развернуто 1 1 2 1 3 2 3 1 1 1 (1 )(1 2 ) ( ) ( ) ( ) (1 )(1 2 ) (1 ( 1) ) 1 , ( ) ( 1) 1 ( ) N N t t t t N t O t N t − − − − − − − = + + + − − − − + + → +В частности, главный член асимптотики имеет вид 1 1 ( ) , 1, ( ) t t t − → + . (Доказательство приведено в Приложении Б. 76 Нахождение дисперсии остаточной наработки. Рассмотрим теперь остаточную дисперсию условной случайной величины ( ) t X T t T t = − . Согласно [133] она выражается следующим образом 2 2 2 2 2 ( ) ( ) ( ) (1 ( )) ( ) ( ) 1 ( ) t t t E X t F x x dx t F В случае распределения Гнеденко-Вейбулла ( , ) GW имеем ( ) ( ) ( ) ( ) 2 2 ( ) 2 ( ) ( ) ( ) ( ) ( ) 2 2 2 ( ) 2 ( ) ( ) ( ) 2 2 ( ) ( ) ( ) ( ) ( ) 2 2 ( t x x z t e e e e dz dx t t x z t z t z e dx e dz t e dz e dx t t x t t t z e e z t dz t t t z z e e zdz t e dz t t + + − − = − = + + + − − = − = − = + Но, ввиду справедливости следующих соотношений ( ) ( ) ( ) t x e dx t t e + − = , ( ) ( ) ( ) t t z z e dz t t e получаем, что 77 ( ) ( ) ( ) 2 2 ( ) 2 ( ) ( ) ( ) ( ) 2 = 2 2 ( ) ( ) t t z t e e zdz t t t t e t z e e zdz t t t t + − = − − = + В силу того, что 1 2 ( ) , ( ) 2 z e zdz t t + − = , получаем формулу для нахождения дисперсии остаточной наработки 1 2 ( ) 2 2 ( ) 2 , ( ) 2 ( ) ( ) 2 t t e t t t t = − − . (44) Асимптотическое представление дисперсии остаточной наработки и коэффициента вариации остаточной наработки Получим асимптотическое представление дисперсии остаточной наработки из предыдущего соотношения. Заметим, что 2 ,( ) ( ) 2 2 2 ,( ) ( ) 2 1 ,( ) t t e t t t = (45) Используя профессиональный пакет символьных вычислений Wolfram Mathematica, получаем следующее разложение 78 2 ,( ) 2 1 2(1 ) (4 11 7 ) 1 1 2 3 4 1 2 3 1 ( ) ( ) ( ) , ( ) t t O t t t t t − − + = + + + + − , верное при 4 1 0 − . В случае деградационных процессов, как мы знаем, 1 . Окончательно, из формулы (45) получаем 1 4(1 ) 2(1 )(3 5 ) 2 1 ( ) 1 2( 1) 2 3 2 2 2 2 ( ) ( ) t O t t t t − − − = + + + − − (46) Изобразим графически полученную зависимость, взяв первые три члена разложения, при значениях параметров масштаба 1 = и формы Рисунок 15. График 2 ( ) t для значений параметров 1, 2 = = . Проверим полученное соотношение на показательном распределении. В этом случае, при значении 1 = 2 2 1 1 ( ) , ( ) t t = = . В силу того, что ( ) 2 , ( ) 2, (1 ) x t t t xe dx t e t + − − = = = + , получаем 2 3 4 5 0.05 0.10 0.15 0.20 0.25 79 1 1 1 1 2 ( ) 2 (1 ) 2 2 2 2 t t t e t t e − = + − − = , что и требовалось показать. На основании полученных соотношений выше найдем асимптотическое представление для коэффициента вариации ( ) C T t V : ( ) 1 1 1 ( ) 1 2 3 2 ( ) ( ) ( ) t C T t O V t t t t − − = = − + + . (47) Графическое представление полученной зависимости, учитывающее первые четыре члена разложения, при значениях параметров масштаба 1 = и формы Рисунок 16. График коэффициента вариации ( ) C T t V для значений 1 = и 1,1 = |