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

УДК 681.518.5

М.Р. Галиаскаров1, И.В. Рудакова2, Л.А. Русинов3

Введение

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

Большинство ПОТП характеризуется большим количеством технологических параметров и плохо описано математическими моделями, особенно для предава-рийных и нештатных режимов. Вследствие этого значительная часть работ, посвященных вопросам мониторинга состояния ПОТП, базируется на статистических методах контроля, к которым, например, относятся многомерные карты Шухарта [1, 2]. Однако, при применении таких карт возникают сложности их построения для задач большой размерности с существенной взаимной корреляцией между переменными.

Для снижения размерности можно использовать линейный метод главных компонент (МГК), позволяющий перейти к новому меньшему ортогональному базису переменных (главных компонент) без существенных потерь в точности статистической модели [3-5]. Кроме того, достоинством МГК является возможность осуществления мониторинга контролем всего двух статистик.

СРАВНИТЕЛЬНЫЙ АНАЛИЗ ХАРАКТЕРИСТИК МОНИТОРИНГА ТЕХНОЛОГИЧЕСКОГО ПРОЦЕССА С ИСПОЛЬЗОВАНИЕМ ЛИНЕЙНОГО И НЕЛИНЕЙНОГО МЕТОДА ГЛАВНЫХ КОМПОНЕНТ

Санкт- П етербургски й государственн ый технологический институт (технический университет) 190013, Санкт-Петербург, Московский пр., д. 26

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

Ключевые слова: Метод главных компонент, мониторинг технологических процессов, керн-метод главных компонент

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

В работе [6] проводилось исследование моделей, полученных на базе линейного и нелинейного МГК, но оценивалась предсказательная способность таких моделей. В данной работе представлены результаты сравнения возможностей использования этих модификаций МГК для мониторинга ПОТП. В качестве объекта мониторинга рассматривается нелинейный технологический процесс пиролиза бензинов.

Нелинейный метод главных компонент

Линейное преобразование, используемое в МГК, имеет вид:

X = ТРт (1)

где X [пхр] - матрица данных процесса из п выборок измерений с p переменными, нормированными к нулевому среднему и единичной дисперсии, Р[рхр] - матрица нагрузок (матрица главных компонент), Т [пхр] - матрица счетов-проекций исходных переменных на главные компоненты [3].

1 Галиаскаров Марат Рамилевич, аспирант каф. автоматизации процессов химической промышленности, e-mail: gm.iamdeamon@gmail.com

2 Рудакова Ирина Викторовна, канд. техн. наук, доцент каф. автоматизации процессов химической промышленности, e-mail: riv-lilu@yandex.ru

3 Русинов Леон Абрамович, д-р техн. наук, профессор, заведующий каф. автоматизации процессов химической промышленности, e-mail: lrusinov@yandex.ru

Дата поступления - 26 сентября 2013 года

Для реализации линейного МГК может быть использовано сингулярное разложение, в частности матрица Х может быть представлена в виде:

X = Ш¥Т (2)

где и содержит левые сингулярные векторы (ц, ¡е[1,р]), V - правые сингулярные векторы, а I - диагональная матрица с элементами ^ - сингулярными значениями (а2 - собственные числа матрицы ХТХ).

Сравнивая выражения (1) и (2), можно видеть, что матрица нагрузок Р в выражении (1) может быть получена из (2), если положить Р = V. Также можно представить матрицу проекций, как Т = и1. Следует отметить, что формула (1) справедлива, если матрица Х удовлетво-р

ряет условию Xх _ д для всех ] = [1...П]. То есть ис-1=1 ""

ходные данные должны быть предварительно центрированы.

Главное достоинство МГК - возможность существенного сокращения размерности модели. Если часть переменных в массиве данных являются коррелированными, то требуется существенно меньшее число главных компонент q (я<р) для объяснения большей части дисперсии данных. При этом часть дисперсии, оставшаяся необъясненной, попадает в подпространство невязок, которое связывают с шумом процесса.

Основная идея нелинейного метода главных компонент (НМГК) состоит в первоначальном отображении входного пространства размерности р в так называемое линеаризующее пространство Р большой размерности И (1"1»р) с помощью нелинейного отображения (Х]еЯ™^Ф(Х])еР"), а затем в вычислении главных компонент уже в пространстве Р [7, 8].

Выборочная ковариационная матрица данных в пространстве Р будет иметь вид: 1 п

Предполагается, что исходные данные централизованы

(3)

(4)

= -Ка (6)

п

где К - матрица скалярных произведений (Грамм-матрица) с элементами к^ = (Ф(х ^)Ф(х к)).

Для центрирования матрицы К используется выражение:

К = К - КЕ - ЕК - ЕКЕ,

(7)

где каждый элемент матрицы Е равен

1

Далее отображение данных в линеаризующем пространстве Р будет эквивалентно нахождению собственных векторов (а1, а2,.. ап) с собственными числа-ми(Л1,>Л2>...> Лп). Размерность задачи можно снизить посредством использования только первых я собственных векторов. При этом коэффициенты^, а2,.., ад) нормализуются так, чтобы (уг, уг) = 1 для г е [1,я].

Тогда проекции t=[tl,...,tq] текущего вектора х на главные компоненты получаются отображением Ф(х) на собственные векторы Уг в пространстве Р:

К = V Ф( х)) = £ а.ф( х, ),Ф( х)} (8)

]=1

Для решения уравнения (6) и при проецировании данных из входного пространства я" в пространство Р можно избежать нелинейного отображения и расчета скалярных произведений в линеаризующем пространстве, вычислением скалярных произведений в исходном пространстве, вводя керн-функции [7-9]:

к(х,у) = (Ф(х), Ф(у)) (9)

Наиболее часто используются следующие керн-функции, удовлетворяющие теореме Мерсера: - полиномиальное ядро

' Д „ ,

X Ф(^) = 0 и Ф - нелинейная функция, которая отоб-к=1

ражает входные векторы из входного пространства я" в пространство Р.

Аналогично линейному МГК главные компоненты находятся, как собственные векторы \еР ковариационной матрицы СР вместе с собственными числами X (вектор V, соотвествующий наибольшему X становится первой главной компонентой)

1 п

XV = С^ = 1X Ф(х г V)) Ф(х J) ^=1

где X - собственные числа (Х>0), ^Ф(хJ,^ - скалярное

произведение Ф(х^ и V. Из (4) следует, что собственные векторы V могут быть представлены линейными комбинациями преобразованных значений входных данных. Так как ^Ф(хj,-у)^ - скаляр, то все значения V с X * 0 лежат

в интервале [Ф(хО,Ф(хп)] и существуют коэффициенты а, такие что

v = X а jФ(х j) (5)

J=l

Умножая на Ф(х¡) слева обе части выражения (4) с учетом выражения (5), получим:

к(х,у) = (х,у) = (хт • у + г)'

- сигмоидное ядро (гиперболический тангенс)

к(х,у) = ишЬ(Рд(х,у) + Р1)

- Гауссово ядро (радиальное базисное ядро)

II II2

х - у

к(х,у) = ехр(

2а2

)

(10)

(11)

(12)

где г, <3, во, Р1, а - должны быть априорно специфицированы пользователем. Полиномиальное и Гауссово ядро всегда удовлетворяют теореме Мерсера, а сигмоидное ядро -только при определенных значениях коэффициентов во и в1.

Для определения количества главных компонент я обычно используется критерий, определяющий долю общей дисперсии исходных данных, объясняемую базисом главных компонент. Иногда считается доля CPV суммы учитываемых моделью МГК собственных чисел X¡ ковариационной матрицы по отношению к сумме всех собственных чисел матрицы СР:

я

XX г

^ г % (13)

СРУ = 100^=1—

п

XX,

1=1

Если эта доля достигает 90 % и более, то считается, что модель адекватно описывает корреляционные связи переменных в процессе.

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

п

дает методов для восстановления данных в первоначальном пространстве. С этим связана основная проблема использования НМГК для мониторинга и диагностики.

Мониторинг состояния процесса на базе НМГК

Мониторинг на базе линейного МГК производится посредством контроля динамики изменения двух статистик: по Т2-статистике контролируются отклонения процесса, объясняемые моделью МГК, а по <-статистике -невязки модели. Каждая из статистик сравнивается с соответствующими порогами, и при их превышении хотя бы одной из статистик принимается гипотеза о появлении нарушения в ходе процесса [4, 5].

Мониторинг на основе НМГК подобен мониторингу на базе линейного МГК в том, что Т2 и <-статистики в линеаризующем пространстве интерпретируются одинаковым образом.

Состояние процесса в пределах модели НМГК контролируется статистикой Т2 - суммой нормализованных квадратичных счетов:

Т2 = ... 1ч]Л-1[1!... Ц]т (14)

где tk рассчитываются по (8), л1 - диагональная матрица из инверсий собственных чисел, соответствующих учитываемым в модели главным компонентам. Пороговые значения Т2-статистики определяются с использованием распределения Фишера:

с = д(п -1) Р

Т Ч,(п—а), а

п — д 4 4

(15)

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

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

Одна из методик вычисления статистики < в линеаризующем пространстве F приведена в [9], где предложено использовать линейный МГК в линеаризующем пространстве Р для вычисления счетов ^ используя все р и только принятые для модели я главных компонент. Тогда для расчета <-статистики используется следующее выражение:

2 р а

о=|фр(х)—фч(х)| =£ 12—2tJ2 (16)

Пороговое значение С< для <-статистики может быть определено посредством приближения взвешенного х2-распределения:

с

о

(17)

где

ё = Ъ/ - весовой параметр, 2 = 2а / - число сте-6 /2a /Ъ

пеней свободы, а и Ь - оцененные среднее и дисперсия <статистики [9]. Использование х2-распределения основано на предположении нормального распределения ошибок, что не всегда верно на практике. Однако, т.к. g и || получаются непосредственно из моментов выборочного распределения данных нормального режима работы процесса, то взвешенное х2-распределение работает хорошо даже в случаях, когда распределение ошибок не соответствует Гауссовому.

Порядок операций при проведении мониторинга

При использовании линейного МГК мониторинг процесса для обнаружения нарушений включает следующие операции.

1. Подготовка системы к работе:

- формирование массива данных Х при нормальном ходе процесса и построение МГК-модели, т.е. определение главных компонент по выражению (2), оценка числа главных компонент я, которые будут учтены в модели, используя (13);

- расчёт матрицы счетов t и по выражениям, аналогичным (15) и (17) - определение пороговых значений Т2 и <- статистик.

2. Проведение мониторинга состоит в:

- вычислении по поступившему новому вектору данных х нового вектора счетов t и значений статистик Т2 и <;

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

В случае НМГК порядок операций при реализации мониторинга аналогичен. Однако на первом шаге, предварительно необходимо выбрать тип используемого ядра, а после формирования массива исходных данных построить Грамм-матрицу и отцентрировать ее, используя выражение (7). Далее найти собственные векторы а из (6) и провести их нормализацию {ак, ак) = 1/Як. По выражению (8) определяются счета t и по (15) и (17) осуществляется расчет пороговых значений Ст и С< .

На втором шаге проведение мониторинга реализуется с учетом необходимости вычисления нового значения к(х,х]) и нового вектора счетов по (9). Новые значения статистик Т2 и < сравниваются с порогами.

Сравнительное исследование методов мониторинга на базе МГК и НМГК

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

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

Б1 - «коксование змеевика радиантной секции», которая характеризуется: ростом давления на входе в радиантную секцию Рвх1 и ростом температуры на выходе из нее Твых1, причем динамика по первому параметру будет опережать отклонение второго из-за маскирующего эффекта системы регулирования;

Б2 - «попадание в змеевик воды с паром разбавления», которую можно выявить по совокупности следующих отклонений: рост давления на стадии закалки Рзиа2, увеличение содержания кислорода Со2 на перетоке печи и

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

Рбвх2-

Эксплуатационные пороги устанавливались по результатам экспертного опроса. Для работы выбраны 2 массива данных, включающих 5 перечисленных выше диагностических параметров, каждый массив содержал 198 временных точек и 15 параметров. Дискретность опроса датчиков составляла 1 минуту.

В качестве ядер НМГК использованы: полиномиальное ядро (10), для которого г = 1 и d = 1, и сигмоид-ное ядро (11), для которого во = 0,005, Р1 = 0,1. Для обеспечения допустимой адекватности модели (90 % -процент объяснения дисперсии данных) достаточно для линейного МГК оставить 7 главных компонент q, а для НМГК - 10, что показано на рисунке 1.

НМГК

[полиномиальное)

Рисунок 1. Выбор числа главных компонент.

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

Для моделирования нештатной ситуации возмущения налагались на Рвх1 и Твых1 с 10 шага, выход переменных за пороговые значения наблюдался на 100 и 107 шагах соответственно. Для моделирования нештатной ситуации Б2 возмущения налагались на Рзиа2, 1=бвх2 и С02 с 50 шага выход переменных за пороговые значения происходил на 170, 150 и 180 шагах соответственно.

Так как исходные данные сильно зашумлены, то массив, на котором выполнялось имитационное моделирование, был подвергнут предварительной фильтрации. Был применен комбинированный фильтр экспоненциального сглаживания и медианы с коэффициентом сглаживания первого фильтра 0,5 и объемом выборки для медианы 10. На рисунке 2 показаны тренды статистик р и Т2 при нормальном протекании процесса. Видно, что статистика Т2 при линейном МГК чревата ложными обнаружениями, для НМГК, та же ситуация наблюдается для р-статистики.

На рисунках 3 и 4 представлены отклики статистик при развитии нештатных ситуаций Б1 и Б2 соответственно. Развитие ситуации начинается с 10-го шага. Числа в рамках соответствуют моментам обнаружения нарушения.

-Св=264г8-

10 20 30 40 50 60 г С =8,17.-т-,-,-_ 70 80 90

кма

10 20 30 40 50 60 70 80 90

: <^=47,78 ]........;.......1.......'

<3 0,02

0,015

0,01

10 20 30 40 50 60 70 80 90 1 " С^-0,02044—|-у-I-Г-Т-Г

лМ!

10 20 30 40 50 60 70 80 90 1

у.

10 20 30 40 50 60 70 80 90 *

10 20 30 40 50 60 70 80 901

Рисунок 2. Тренды статистик Q и Т при нормальном ходе процесса пиролиза: а - МГК, б - НМГК (Полиномиальное ядро), в - НМГК (Сигмо-идное ядро).

32,92

!

'Й

И 20 30 40 50 60 70 80 90

С^-8,17

С, г 41 73

20 30 40 50 60 70 80 901

40 50 60 70 80 90

О 0,025 0,02 0,015 0,01 0.005

<^=0,02044

20 30 40 50 60 70 80 901

40 50 60 70 80 90 1

Рисунок 3. Тренды статистик <3 и при развитии первой нештатной ситуации: а - МГК, б - НМГК (Полиномиальное ядро), в - НМГК (Сигмо-идное ядро).

Г,-264,8 *-

* :

10 20 50 40 й 60 70 8 90 100

С, 47.78

10 20 30 40 50 60 70 80 ¿,90

<>-47.78

У !

|

10 20 30 40

60 70 80 90

10 20 30 40 50 60 70

Рисунок 4. Тренды статистик (? и 1*при развитии второй нештатной ситуации: а - МГК, б - НМГК (Полиномиальное ядро), в - НМГК (Сигмо-идное ядро).

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

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

Литература

1. Царев Ю.В. Статистические методы контроля и управления качеством. Контрольные карты: учебно-методическое пособие. Иваново: ИГХТУ, 2006. 160 с.

2. Bersimis S. Psarakis S., Panaretos J. Multivariate Statistical Process Control Charts: An Overview. // Qual. Re-liab. Engng. Int. 2007. V. 23. P. 517-543.

3. Эсбенсен К Анализ многомерных данных. Черноголовка: Изд-во ГЕОХИ РАН, 2005. 158 с.

4. Bishop C.M. A review of process fault detection and diagnosis. N.Y.: Prentice Hall, 2000. 175 p.

5. Venkatasubramanian V, Rengaswamy R., Yin K, [et al] A review of process fault detection and diagnosis. Quantitative model-based methods // Computers and Chemical Engng. 2003. V. 27. P. 293-346.

6. Caoa L.J., Chuab KS., Chongc W.K [et al] A comparison of PCA, KPCA and ICA for dimensionality reduction in support vector machine // Neurocomputing. 2003. V. 55. P.321-336.

7. Schölkopf B., Smola A.J., Müller K.-R. Nonlinear component analysis as a kernel eigenvalue problem. // Neural Computation. 1998. V. 10. P.1299-1319.

8. Choi, S.W., Leeb Ch., Leeb J-M. [etal] Fault detection and identification of nonlinear processes based on kernel PCA // Chemometrics and Intelligent Laboratory Systems. 2005. V. 75. P. 55-67.

9. Lee J. Yoo Ch. Y, Choi S.W. [et al]. Nonlinear process monitoring using kernel principal component analysis. // Chemical Engng Sci,2004. V. 59. P.223-234.