МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ, ЧИСЛЕННЫЕ МЕТОДЫ И КОМПЛЕКСЫ ПРОГРАММ
УДК 621.452:519.217
Г. Г. Куликов, В. Ю. Арьков, А. И. Абдулнагимов
НЕЧЕТКИЕ ИЕРАРХИЧЕСКИЕ МАРКОВСКИЕ МОДЕЛИ ДЛЯ КОНТРОЛЯ И ДИАГНОСТИКИ ГАЗОТУРБИННЫХ ДВИГАТЕЛЕЙ
Предложена методика определения состояния газотурбинных двигателей на основе нечетких иерархических марковских моделей. Обсуждается использование нечетких систем в качестве инструмента интерполяции матрицы вероятностей переходов между базовыми точками. Показана возможность использования формул Байеса как инструмента нечеткой логики при идентификации и моделировании марковских моделей. Газотурбинные двигатели; система автоматического управления, контроля и диагностики; системная безопасность; идентификация систем; марковские модели; матрица вероятностей переходов; нечеткая логика; анализ иерархий
ВВЕДЕНИЕ
Безопасность представляет собой критическую проблему летательных аппаратов без решения которой их эксплуатация недопустима. Согласно статистике, в последние годы большинство происшествий в авиации связаны с поздним обнаружением отказов в системах самолета. В связи с этим с каждым годом повышаются требования к обеспечению безопасности полета, которые требуют разработку новых методов и алгоритмов управления, контроля и диагностики за сложными объектами.
Безопасность полета самолета зависит в значительной степени от надежной работы силовой установки, частью которой является бортовая система управления и контроля типа FADEC. Силовая установка должна сохранять работоспособность даже при наличии отказов; для этого «опасный» отказ переводится в «безопасный» за счет реконфигурации системы. Примером является снижение режима работы двигателя при превышении допустимых значений параметров, вывод из помпажа и др.
Одним из перспективных направлений в обеспечении системной безопасности является разработка интеллектуальных информационных систем с алгоритмами контроля и диагностики на основе нечетких иерархических марковских моделей.
Для описания сложных динамических систем с учетом случайных возмущений предлагается использовать модели Маркова. С точки зрения динамического моделирования, марков-
Контактная информация: 8(347) 273-78-23 Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (гранты №12-08-31279, 12-08-97027)
ские цепи представляют собой непараметрическую модель (многомерную гистограмму) в форме многомерного распределения, которое, в общем случае, может изменяться во времени. Управляемые цепи Маркова учитывают влияние управляющих входных сигналов на поведение динамической системы [1]. В то же время остается открытым вопрос о переходе от стохастических моделей в форме управляемых марковских цепей к детерминированным динамическим моделям (таким, как спектры и весовые функции), являющимся общеупотребительным экспертным инструментом анализа и синтеза систем автоматического управления, контроля и диагностики (САУКиД) [2].
Марковские цепи позволяют естественным образом объединить в модели как динамические свойства объекта (спектральные и весовые характеристики), так и параметры распределения наложенных возмущений. Кроме того, в этом случае появляются возможность «извлекать» из стохастической модели информацию о свойствах моделируемой системы в желаемой форме.
1. ИЕРАРХИЧЕСКАЯ МОДЕЛЬ СОСТОЯНИЯ СИЛОВОЙ УСТАНОВКИ
Для разработки интеллектуальной системы управления необходимо формализованное представление силовой установки в виде иерархической структуры как сложной системы, состоящей из подсистем и элементов со встроенными функциями контроля (ВФК).
Для этого проводится декомпозиция силовой установки на независимые подсистемы различного уровня иерархии по структурным и функциональным признакам, например:
• САУКиД (FADEC);
• гидромеханическая исполнительная часть САУКиД;
• топливные (основную и форсажную) камеры;
• система запуска;
• масляная система;
• дренажная система и др.
На рис. 1 приведена обобщенная схема иерархической структуры силовой установки в соответствии с требованием 1. Предполагается, что в каждом элементе или подсистеме ГТД заложена ВФК физически или на логическом уровне.
При построении САУКиД проводится классификация, формализация и представление процессов контроля состояния и диагностики отказов для основных частей двигателя (компрессор, камера сгорания, турбины, сопло, и т. д.) и его систем (управления, контроля, подачи топлива и т. д.) [3]. Эти процессы можно представить в форме комплексных дискретных марковских цепей (диагностических марковских таблиц), которые позволяют провести анализ состояния силовой установки.
Марковская матрица вероятностей переходов для моделирования отказов и их последствий имеет универсальную структуру для всех трех уровней декомпозиции системы:
• система в целом (силовая установка);
• конструктивные блоки;
• элементы.
Иерархическая структура предложенной марковской модели строится в обобщенном пространстве состояний, где двоичные флаги отказов и реконфигурации образуют слово состояния элемента, блока, силовой установки, что позволяет проследить развитие отказа и процесс его диагностики. В ходе проектирования системы контроля, главным образом рассматривается область единичных отказов. Предложенная марковская модель дает возможность представить систему с множественными отказами и последовательностями отказов. Верхний уровень декомпозиции слова состояния силовой установки отражает все возможные отказы силовой установки в целом. Принимая во внимание, что моделирование всех множественных отказов и последовательностей отказов потребовало бы слишком много времени, выбираются только наиболее вероятные комбинации на основе статистики эксплуатации и экспертных знаний (рис. 2).
Марковская модель представляет собой вероятностную модель, отражающую семейство процессов диагностики [3]. Каждый процесс исследуется на полунатурном испытательном стенде [4, 5], где реальное оборудование системы управления и контроля FADEC работает в замкнутом контуре с имитационной моделью двигателя. Поскольку для заполнения марковской матрицы необходимо провести большое количество однотипных испытаний САУКиД с моделированием различных отказов и их комбинаций, предлагается проводить такую программу тестов и заполнения матрицы в автоматизированном режиме.
Рис. 1. Декомпозиция силовой установки
Шаг Т,:
Г
-Л.
Шаг Т+
-Л.
Единичные
Множественные отказы
Для оценки влияния каждого отказа подсистемы ГТД на всю силовую установку летательного аппарата используется нечеткий алгоритм идентификации состояния системы и метод анализа иерархий [6]. Для каждой подсистемы и элемента ГТД строятся матрицы суждений на основе экспертных знаний и статистики, где вычисляются коэффициенты влияния отказов.
В рамках разрабатываемой методологии рассматриваются не только одиночные отказы (как принято в теории и практике проектирования систем управления и контроля), но и комбинации отказов, а также последовательности (цепочки) так называемых «следственных» отказов. К тому же, с помощью марковской модели можно выявлять и контролировать «ложные отказы».
На базе нечетких правил определяются состояния элементов или подсистемы на каждом уровне иерархии. Отказные состояния классифицируются в зависимости от степени деградации на «незначительное», «сложное», «опасно сложное», «катастрофическое».
Состояние всей силовой установки определяется суммой произведений коэффициентов влияния отказов на степень деградации элемен-
та или подсистемы, после чего строится вероятностная модель для всей системы.
В соответствии с формализованным представлением, можно определять динамическую модель состояния всей силовой установки как снизу вверх, так и сверху вниз. Данная модель состояний позволяет оценивать «расстояние» до критической ситуации. Визуальный анализ тренда дает возможность оценить резерв времени до наступления критической ситуации, и, таким образом, позволяет заблаговременно планировать действия экипажа.
2. НЕЧЕТКИЕ ИЕРАРХИЧЕСКИЕ МАРКОВСКИЕ МОДЕЛИ
Пусть ^2, ... , ¥п - совокупность отказов,
которые встречаются в подсистемах ГТД некоторого уровня иерархии.
Необходимо определить веса шь ш2, ... , юп и их влияние на системы ГТД следующего уровня и силовой установки в целом. Весовые коэффициенты отражают степень влияния различных отказов на безопасность полета. В частности, эксперт может оценить три вида отказов, например:
1. Обрыв первой или второй обмотки датчика вращения вала;
2. Отказ первого канала системы управления;
3. Отказ резервного канала системы управления.
В первом случае - рассмотрим датчик вращения вала с двумя идентичными индуктивными обмотками. После обрыва первой обмотки система управления обнаружит, что электрическая цепь разомкнута, и реконфигурирует систему так, что она будет использовать вторую обмотку для измерения параметров.
Во втором случае, когда оборваны обе обмотки датчика, это говорит о том, что датчик не может быть использован далее для измерений параметров. Система управления обнаруживает данный отказ и переводит систему на другой контур регулирования. Вместо измерения скорости вращения вала система управления измеряет давление в компрессоре, при этом управляя тем же расходом топлива.
В третьем случае, при отказе первого канала системы управления, что определяется встроенным тестовым алгоритмом, система переходит на второй (резервный) канал.
Когда отказывают оба канала цифровой системы управления, управление осуществляется гидромеханической системой.
Эти примеры показывают влияние каждого отказа на состояние всей силовой установки. Поэтому эксперт в состоянии оценить относительную значимость вышеперечисленных отказов. Такая классификация отказов делает возможность проведения формального описания состояния системы методом анализа иерархий.
Количественные суждения о значимости отказов строятся для каждой пары отказов (Fu Fj) и представляются матрицей A размера «х«.
A=(alj), (i, j = 1,2,3), где число aij показывает значимость отказа Fi по сравнению с Fj.
Элементы aij определены по правилам:
1. Если aij = а, то aji = 1/а, а Ф 0.
2. Если суждения таковы, что отказ Fi имеет одинаковую с Fj относительную важность, то aij = 1, aji = 1, в частности aii = 1 для всех i.
Таким образом, получаем обратно-симметричную матрицу A [6]
A =
1/ а.
!/а1п 1/а2п
1
После представления количественных суждений о парах F]) в числовом выражении
через числа Оу задача сводится к тому, чтобы п возможным отказам Fl, F2,..., Fn поставить в соответствие множество коэффициентов влия-яния отказов юь ю2,..., юп, которые соответствовали бы суждениям о значимости отказов систем ГТД.
Если суждение эксперта совершено при всех сравнениях, то матрицу А называют согласованной.
Если А - матрица значений парных сравнений отказов, то для нахождения вектора приоритетов для классификации отказов необходимо найти вектор ю, который удовлетворял бы критерию
Аю Хтахю.
Имеем ю - собственный вектор матрицы А и ^тах - максимальное собственное значение, близкое к порядку матрицы п, тогда
F
A = 2
F
Ю1/ 'Ю1 Ю1/ 'Ю
Ю2 /ю1 Ю2 /ю
Юп /ю1 Юп /ю
®l/ Юп
Ю2/ Юп
®п/Юп
= 1
Заметим, что малые изменения в аі}- вызывают малое изменение Атах, отклонение последнего от п является мерой согласованности. Оно позволяет оценить близость полученной шкалы к основной шкале отношений. Следовательно, индекс согласованности
(^тах - П)/(П - 1)
рассматривается как показатель «близости к согласованности». В общем случае, если это число < 0,1, то можно быть удовлетворенным суждениями о значимости отказов.
Для того чтобы использовать марковские модели силовой установки в качестве динамической системы, рассматривается взаимосвязь матрицы вероятностей переходов с основными статистическими свойствами сигналов. В каче-
1
а
а
1п
1
а
F
F
2
стве исходной математической модели объекта управления рассматривается описание в пространстве состояний в виде стохастического разностного уравнения
ХЦ+ \)=АХЦ)+Щг), (1)
где X е В? - 5 -мерный вектор состояний; Л и /■ -матрицы размерности (п У п) и (п У г) соответственно; е К - вектор независимых случайных переменных. Уравнение (1) по определению является процессом Маркова, поскольку будущее состояние системы Х{!+1) определяется предшествующим состоянием Х{!) и не зависит от предыстории процесса [7, 8].
Введение дискретизации по уровню позволяет перейти от марковского процесса к цепи Маркова. При условии стационарности процесса £,(0 цепь будет однородной. Такая цепь описывается с помощью стохастической матрицы вероятностей переходов Р размерностью {тхт), где т - количество состояний цепи. Каждый элемент матрицы /'(; представляет собой вероятность перехода системы из состояния в состояние Х}- за время АТ:
Р = РгоЬ{Х (?) = X,., Х{г +1) = Х]} УneN;
___' (2)
Ру =1, 1=1,т
1
Условие (2) означает, что матрица/* должна быть стохастической и определять полную систему событий. Сумма элементов в каждой строке матрицы должна быть равна единице. В случае, когда доступны наблюдения не только состояния X, но и управления £/, модель (1) принимает вид:
ХЦ + 1) = АХ{ 0 + ВЩ 0 + Щ 0- (3)
При анализе стационарных случайных процессов текущие вероятности состояний со временем приближаются к стационарным (финальным) значениям. Поэтому оценка плотности распределения в виде гистограммы должна совпадать с оценкой распределения как вектора финальных стационарных вероятностей.
Для марковского моделирования размер интервала дискретизации Ахсущественно больше шага квантования по уровню при оцифровке. В качестве события принимаем принадлежность значения х(?п) 7-му интервалу дискретизации:
х(0 = Хг : х(0е [Хг -Ах/2,Хг + Ах/2].
Стохастическая динамика объекта, описываемого уравнением (3), соответствует управляемой цепи Маркова первого порядка.
Рассмотренные свойства относятся к цепи Маркова первого порядка (3) и могут быть обобщены на случай управляемой марковской цепи высокого порядка. Системы более высокого порядка могут быть описаны сложными цепями Маркова либо преобразованы к простым цепям путем введения векторного описания состояний.
Марковские модели стохастической динамики сложных систем можно строить по экспериментальным данным. Размер матрицы Р определяется априорной информацией о порядке модели объекта (1) и числом интервалов дискретизации Ах и А и. Вероятности переходов оцениваются как относительные частоты соответствующих событий.
Статистическая оценка переходных вероятностей управляемой марковской цепи сводится к подсчету соответствующих событий за время наблюдения и последующему вычислению элементов матрицы Р по формуле:
2Х
У=1
где Nijk - число событий вида {Х(?„) = Хь X{t„vО =Xj, U{tn) = Uk}, а знаменатель соответствует числу событий вида {Х(?„) = X, U{tn) = Uk}. Таким образом, при любой комбинации состояния Xj и управления Uk получаем полную систему событий, заключающихся в переходах в состояния Xj. Нормализация (4) делает матрицу Р стохастической. В результате набор вероятностей по каждой строке Ру описывает полную систему событий, для которых сумма вероятностей равна единице:
т
2Л=‘-
У=1
Идентификация марковской цепи аналогична построению гистограммы. Гистограмма представляет собой результат непараметрического оценивания функции плотности распределения в виде вектора вероятностей состояний X.
Распределение вероятностей нового состояния системы X = {Prob(X)} при заданных^ и Uk определяется соответствующей строкой матрицы/’ (рис. 3):
¿ = P\lk = {P1]kJ = l-m,i =
(5)
= const,к = const} = {рлк ...Р1тк}.
Моделирование координат состояния по марковским моделям производится с использо-
ванием метода Монте-Карло и метода функционального преобразования. Для этого на каждом шаге моделирования 4 определяется текущее состояние X, и текущее управление Затем находят распределение (5) и вычисляют интегральную функцию распределения Г:
л
Г (х) = | р( x)dx,
Г (х1) = РгсЬ{ X (/„ ) < X; } = ^Рм.
(6)
Р; = Рг сЬ{X, ® X, }
1. Исправное состояние
Рз1
г
Рзз
Рис. 3. Матрица вероятностей перехода силовой установки в течение полета
Далее, с помощью генератора случайных чисел У с равномерным распределением и функционального преобразования Г"1 формируется случайная величина X(/n+1) с заданным законом распределения (6):
X(tn+l) = ^(У(4)).
Как показывает анализ, динамическая марковская модель адекватно передает не только форму распределения, но и стационарные (корреляционные) характеристики выходного сигнала [9].
3. КОСВЕННАЯ ОПТИМИЗАЦИЯ ОЦЕНОК МАРКОВСКИХ МОДЕЛЕЙ
При оценке марковских моделей из собранных данных (работы силовой установки), возникает проблема оптимального выбора разрешающей способности. На рис. 4 видно, что при высокой разрешающей способности возрастает случайная погрешность - дисперсия. В этой области график становится зашумленным. При низком разрешении производится усреднение, приводящее к росту систематической ошибки, -
смещения. Здесь график становится гладким, но оценка существенно отклоняется от своего точного значения.
В работе предлагается использовать поиск оптимального разрешения по косвенному графическому критерию. Оцениваются формы распределения с различным разрешением и степенью перекрытия интервалов группирования.
Такое представление матрицы вероятностей переходов и выбор оптимальной длины интервала группировки и перекрытия интервалов гистограммы позволяют найти компромисс между случайными и систематическими ошибками.
Предварительный теоретический анализ показывает, что искомая оптимальная величина разрешающей способности существует. Однако измерить величину суммарной погрешности практически невозможно, поскольку «истинная» модель объекта неизвестна [10]. Поставленная задача может быть решена только косвенным путем.
Ах
Г^тгЧ
Рис. 4. Компромисс между случайными и систематическими ошибками
Не имея истинных оценок, при использовании многократного оценивания с различными интервалами группировки перекрытия интервалов, можно определить область компромисса по графику.
Марковская модель как стохастическая система может передавать основной вид нелинейности и структуру исследуемой модели [7, 8]. Размерность марковской модели определяется структурой детерминированной модели и видом возможной нелинейности. На рис. 5 приведен пример расхода топлива на установившемся режиме работы ГТД. Видно, что расход топлива меняется со временем в форме линейного тренда. Колебания намного больше уровня квантования в цифровой системе (маленькие колебания в средней области графика).
Время
Рис. 5. Реализация исходного процесса (От) расхода топлива
Такое поведение расхода топлива объясняется развивающейся нелинейностью
в гидромеханической части системы. Возможный тип нелинейности в этом случае - сухое трение или зона нечувствительности в гидроцилиндре, который приводится в движение гидравлическим давлением. В этом случае цилиндр теряет свою чувствительность к небольшим изменениям при гидравлическом давлении, и он может только перемещаться под более высоким давлением. Однако такая статическая ошибка компенсируется обратной связью в системе управления со случайными возмущениями.
В результате можно наблюдать отсутствие статической ошибки и определенного типа динамической ошибки из-за развивающегося отказа в гидромеханической системе. На данном этапе рассматривают систему, работающую в обычном режиме, тогда как предложенный подход марковского моделирования позволяет обнаруживать развивающиеся отказы.
Построение матриц вероятностей переходов состояний для САУКиД ГТД проводится в среде технических вычислений МаАаЬ.
Моделирование марковских процессов включает в себя следующие шаги: дискретизацию состояния и управления по уровню, статистическое оценивание матрицы вероятностей переходов, задание начального распределения вероятностей состояний, вычисление функции распределения по данному управлению и пошаговую генерацию состояний при заданном распределении.
На рис. 6 показано трехмерное представление матрицы вероятностей переходов из реализации исходного процесса. Можно уви-
деть, что главный пик соответствует медленной основной тенденции временного ряда, а появление нескольких соседних пиков меньшей амплитуды может быть вызвано развивающейся нелинейностью.
На практике предлагается использовать параллельные вычисления из-за высокой вычислительной нагрузки во время оценок матриц вероятностей переходов [10].
Используется принцип распараллеливания по данным, поскольку вычисление элемента матрицы переходных вероятностей не зависит от соседних элементов. В работе применяется параллелизм на уровне подсчета строк, столбцов и блоков матрицы, где различные исполняемые модули выполняют свой программный код над своим набором данных.
Рис. 6. Трехмерное изображение сечения матрицы вероятностей переходов
Данный подход позволяет сбалансировать загрузку вычислительных модулей кластера или суперкомпьютера.
ЗАКЛЮЧЕНИЕ
Рассматривается технология марковского моделирования для обнаружения отказов, основанная на выявлении нелинейной динамики и дальнейшем анализе иерархического состояния газотурбинной силовой установки. При переходе к более высокому уровню иерархии, данные нижних уровней агрегируются, чтобы получить обобщенные показатели состояния системы, подходящие для индикации в кабине пилота и поддержки принятия решений со стороны пилота или экипажа. С этой целью метод анализа иерархий был объединен с марковским представлением и использованием «слов состояний».
СПИСОК ЛИТЕРАТУРЫ
1. Куликов Г. Г., Арьков В. Ю., Брей-
кин Т. В. Марковское моделирование динамических объектов для полунатурных испытаний // Известия РАН. Теория и системы управления. 2000. № 2. С. 124-128.
2. Isermann R. Fault-Diagnosis Systems. Berlin, New York, Springer-Verlag, 2006. 475 р.
3. Dynamic modelling of gas turbines: identification, simulation, condition monitoring, and optimal control / G. Kulikov G. [et al.]; G. Kulikov and H. Thompson, eds., London, New York, Springer, 2004.
4. Breikin T. V., Arkov V. Y., Kulikov G. G. Application of Markov chains to identification of gas turbine engine dynamic models // International Journal of Systems Science. February 2006. Vol. 37, No. 3. P. 197-205.
5. Куликов Г. Г., Арьков В. Ю., Абдулнагимов А. И. Технология полунатурных испытаний интегрированных систем управления и контроля авиационных ГТД // Известия вузов. Авиационная техника. 2008. № 1. С. 37-40.
6. Методология полунатурного комплексного функционального моделирования ГТД и его систем / Г. Г. Куликов [и др.] // Вестник УГАТУ. 2009. T. 13, № 2 (35). C. 88-95.
7. Куликов Г. Г., Арьков В. Ю., Абдулнаги-мов А. И. Нечеткие марковские модели систем автоматического управления и контроля // Интеллектуальные системы управления / под ред. акад. РАН С. Н. Васильева. М.: Машиностроение, 2010. С. 154163.
8. Saaty T. L. Analytic Hierarchy Process: Planning, Priority Setting, Resource Allocation, McGraw-Hill, New York, London. 1980.
turbine power plant // Proc. IFAC Conf. on Control Methodologies and Technology for Energy Efficiency, Portugal, 2010.
10. Куликов Г. Г., Арьков В. Ю., Абдулнагимов А. И. Идентификация марковских моделей нестационарных динамических объектов на основе параллельных вычислений // Параллельные вычисления и задачи управления: тр. 5-й Международн. конф. PACO’2010. М.: Ин-т проблем управления им. В. А. Трапезникова РАН, 2010. С. 343-350.
ОБ АВТОРАХ
Куликов Геннадий Григорьевич, проф., зав. каф. автоматизированных систем управления. Дипл. инж. по автоматиз. машиностроения (УАИ, 1971). Д-р техн. наук по системн. анализу, автоматическ. управлению и тепловым двигателям (УАИ, 1989). Иссл. в обл. АСУ и управления силовыми установками летательн. аппараторв.
Арьков Валентин Юльевич, проф. той же каф. Дипл. инженер по промышл. электронике (УАИ, 1986). Д-р техн. наук по системн. анализу и управлению (УГАТУ, 2001). Иссл. в обл. спектральн. анализа и идентификации динамическ. моделей авиационных двигателей, систем искусственного интеллекта.
Абдулнагимов Ансаф Ирекович, ст. преп. той же каф. Дипл. магистр техники и технологии (УГАТУ, 2007). Канд. техн. наук по системн. анализу и управлению (УГАТУ, 2012). Иссл. в обл. автоматическ. управления, идентификации и системной безопасности авиационных двигателей.
9. Kulikov G., Arkov V., Abdulnagimov A.
Markov modelling for energy efficient control of gas