Идентификация промышленных объектов по данным их нормального функционирования (часть 1) Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 621.314

ИДЕНТИФИКАЦИЯ ПРОМЫШЛЕННЫХ ОБЪЕКТОВ ПО ДАННЫМ ИХ НОРМАЛЬНОГО ФУНКЦИОНИРОВАНИЯ (ЧАСТЬ 1)

© А.В. Баев1, А.А. Мельник2

Иркутский государственный технический университет, 664074, Россия, г. Иркутск, ул. Лермонтова, 83.

Рассмотрена методика идентификации объектов и систем на математических моделях с помощью пакета Sistem Identifcation Toolbox матричной лаборатории MatLab + Simulink. В результате проведенного эксперимента определены оптимальные статические и динамические характеристики промышленных объектов по данным их нормального функционирования. Предложены рекомендации по применению разработанной методики для идентификации промышленных объектов и систем. Ил. 11. Библиогр. 3 назв.

Ключевые слова: идентификация; система управления; оптимизация; математическая модель; эксперимент; статические и динамические характеристики.

IDENTIFICATION OF INDUSTRIAL FACILITIES BY THEIR NORMAL OPERATION DATA (PART 1) A.V. Baev, A.A. Melnik

Irkutsk State Technical University, 83 Lermontov St., Irkutsk, 664074, Russia.

The article deals with the method of object and system identification on mathematical models using the System Identification Toolbox package of MatLab + Simulink matrix laboratory. Conducted experiment results in the determination of optimal static and dynamic characteristics of industrial facilities by the data of their normal operation. Recommendations on the application of the developed methodology for industrial facility and system identification are proposed. 11 figures. 3 sources.

Key words: identification; control system; optimization; mathematical model; experiment; static and dynamic characteristics.

В настоящее время большинство объектов управления оснащено микропроцессорной техникой, которая позволяет получать информацию о состоянии объекта, обрабатывать ее по определенным алгоритмам и хранить довольно продолжительное время. Эта информация помогает проводить анализ хода технологического процесса, своевременно выявлять нарушения в технологии и причины их появления, при необходимости вносить корректировку в управление. Наряду с этим возможно использование такой информации для определения динамических характеристик как объекта, так и системы регулирования, т.е. проведения идентификации. Попытаемся ответить на вопрос об идентификации с точки зрения практики, используя стандартные средства, в частности пакет System Identification Toolbox матричной лаборатории MatLab+Simulink. Во-первых, всегда ли возможна идентификация, во-вторых, если возможна, то при каких условиях.

Идентифицируемый объект (систему) принято представлять в виде, представленном на рис. 1, где t - время; u(t) - входной сигнал; - теоретический выход объекта; у(^ - наблюдаемый выход объекта; e(t) - эквивалентная случайная помеха, отражающая действие шумов на объект.

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

Можно ли в этих условиях по имеющейся информации в виде временных рядов входного сигнала и

► y(t)

Рис. 1. Схема определения динамических характеристик статистическими методами

1Баев Анатолий Васильевич, кандидат технических наук, профессор кафедры автоматизации производственных процессов, тел.: (3952) 405243, e-mail: [email protected]

Baev Anatoly, Candidate of technical sciences, Professor of the Department of Automation of Production Processes, tel.: (3952) 405243, e-mail: [email protected]

2Мельник Альбина Аполлинарьевна, аспирант, тел.: (3952) 405878, e-mail: [email protected] Melnik Albina, Postgraduate, tel.: (3952) 405878, e-mail: [email protected]

реакции системы на этот сигнал определить динамические характеристики объекта? Рассмотрим упрощенную структуру одноконтурной системы автоматического регулирования (рис. 2).

Согласно схеме, представленной на рис. 2, на вход 1 поступает задающее воздействие, которое будет содержать случайную составляющую только в случае следящей системы или в системе соподчиненного управления (в системе стабилизации случайная составляющая отсутствует). По входу 2 вероятность появления случайной составляющей наибольшая, так как он представляет собой поток регулирующей среды, протекающей через регулирующий орган. В этом случае любые изменения параметров регулирующей среды, которые носят случайный характер, оказывают влияние на переходные процессы в системе. Вход 3 до некоторой степени является условным: сюда методом структурных преобразований перенесены воздействия по всем остальным каналам объекта, включая шум измерительных устройств.

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

По каналу выход у(1) - х ^(?) (1 канал воздействия)

W3c (pУ

Wp - р (Р) (Р)

об

1+WP - р (Р)• ( р)'

об

По каналу выход y(t) - f 1(t) (2 канал воздействия)

W3C (Р У

Wo6(Р)

1+Wp - р (Р) • (Р)

об

По каналу выход y(t) - f 2(t) (3 канал воздействия)

W3C (Р) = 1

1+Wp- р (p) • Wo6 (p)

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

Wзс(р) при условии точного определения переда-

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

Wo6(Р) =

Wo6(Р) =

1

W3C (Р)__

I-W3C(Р) Wp-p(p)'

W3C (Р) 1- W3C (Р) Wр-р (p)'

Wo6(Р) = '

1-W3C(P) 1

W3C (Р) Wp-p (ру

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

Можно ли получить передаточную функцию объекта регулирования, если входом объекта полагать

сумму f1(t) + u (t) ? По [1] эта операция считается

ошибочной, но автор рассматривает операцию параметрической идентификации с точки зрения построения адаптивной системы. Не вдаваясь в теоретические выкладки, рассмотрим вопрос идентификации с помощью пакета System Identification Toolbox в среде MatLab+Simulink [2, 3] в идеальных условиях с помощью модели системы автоматического регулирования. С этой целью в Simulink^ построим модель, позволяющую менять воздействие по различным каналам (рис. 3).

В этой схеме подпрограммой Subsystem определен объект регулирования, а Subsysteml - генератор случайного сигнала. Блоки Gain3, Gain4 и Gain5 предназначены для включения или отключения отдельных частей схемы с целью придания определенных свойств системе, блоки Const и Constl - для изменения задания или внесения возмущения на систему.

2 /2(0 3

О 1 г

¿У УУ

y(t)

f 1(t)

Рис. 2. Структурная схема одноконтурной системы регулирования

Рис. 3. Упрощенная обобщенная структурная схема модели системы для идентификации объекта управления с входным сигналом по различным каналам

Рис. 4. Реакция объекта на скачок и случайное воздействие

В случае, когда значения Const1 и Gain3 равны 0, исследуется объект Subsystem с фиксацией входного сигнала, представляющего собой сумму постоянного воздействия и случайной составляющей, и выходного сигнала с блоками передачи данных Out1 и Out2 соответственно. При других сочетаниях можно рассмотреть замкнутую систему автоматического регулирования.

Для анализа получаемых данных предварительно снимаем динамическую характеристику объекта 2 при воздействии типа скачка 1 без случайной составляющей (естественно, это практически невозможно сделать на реальном объекте) (рис. 4).

Объект устойчивый, обладает незначительным запаздыванием т (приблизительно 5 с) и постоянной времени Т около 55 с. При наложении на скачкообразное воздействие случайной составляющей 3 полу-

чаем реакцию объекта 4.

Для получения математической модели объекта воспользуемся пакетом расширения System Identification Toolbox системы MatLab+Simulink. Предварительную обработку данных проведем в среде GUI пакета. Порядок работы подробно описан в [2].

Для получения модели объекта проведем следующие эксперименты в соответствии с рис 3:

- на вход объекта подается только случайный сигнал от генератора (const=0, const1=0, Gain3=0, Gain4=0, Gain5 =5);

- на вход объекта подается постоянная составляющая и случайный сигнал от генератора (const=10, const1=0, Gain3=0, Gain4=0, Gain5 =5);

- на вход объекта подается постоянная составляющая, случайный сигнал от генератора и случайный процесс замкнутой системы (const=10, const1=0, Gain3=1, Gain4=0, Gain5 =5).

Покажем подробно порядок обработки данных для первого случая. В командном окне определим массивы данных u1 и y1, соответственно входной и выходной сигналы объекта, регистрируемые с помощью блоков передачи данных Out1 и Out2. Командой ident входим в среду GUI пакета идентификации - открывается окно ident:Untitled (рис. 5).

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

Так как структура модели заранее неизвестна, то в окне Parametric Models (а) выбираем структуру State Space и в окне Model Order Selection (б) определяем порядок модели (рис. 7).

Рис. 5. Окно графического интерфейса System Identification Toolbox с результатами построения моделей

И Time Plot: ul->yl_[ ■=■ I ^ [■¿^■J

File Options Style Channel Help

Input and output signals

_i_i_i_i_

0 100 200 300 400 500 Time

Рис. 6. Реакция объекта y1 на ограниченный по модулю входной случайный сигнал u1

В этом же окне Parametric Models можно выбирать следующие структуры: ARX, ARMAX, OE, BJ. Остановимся на методе Бокса-Дженкинса (BJ) как наиболее универсальном. Задаваясь структурой модели в окне Parametric Models, с помощью команды Estimate отправляем в окно идентификации, где в окошечках указывается ее структура. Включив Model output, в левом окне находим переходные процессы в соответствии с выбранной структурой, а в правом окне - степень совпадения переходного процесса с исходным (в нашем случае виден только один переходный процесс, так как степень совпадения приближается к 100%) (рис. 8).

Таким образом, в окне Data/model Info получена модель вида

y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t), где B(q) = 0.002026 qA-4 + 0.007507 qA-5 + 0.001691 qA-6;

C(q) = 1 + 0.09713 qA-1 + 0.001784 qA-2; D(q) = 1 - 2.692 qA-1 + 2.392 qA-2 - 0.6996 qA-3 +1.122e-007 qA-4;

F(q) = 1 - 1.697 qA-1 + 0.7024 qA-2 - 3.934e-009 qA-3 + 1.518e-007 qA-4.

Для сравнения проведем параметрическую оценку по методике Бокса-Дженкинса из командного окна системы MatLab, используя те же данные входного и выходного сигналов: >> z=[y1 u1];

>> th=bj(z, [3 2 4 4 4]); present(th); Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)

B(q) = 0.002026 (±1.103e-012) qA-4 + 0.007507 (±1.396e-011) qA-5

+ 0.001691 (±4.897e-011) qA-6 C(q) = 1 - 0.00327 (±0.0004342) qA-1 + 0.007925 (±0.0001484) qA-2

D(q) = 1 - 1.993 (±0.000429) qA-1 + 1.165 (±0.0007921) qA-2

- 0.1692 (±0.0003677) qA-3 + 1.506e-005 (±9.415e-007) qA-4

F(q) = 1 - 1.697 (±6.888e-009) qA-1 + 0.7024 (±1.312e-008) qA-2

- 5.283e-010 (±7.36e-009) qA-3 - 5.728e-010 (±1.106e-009) qA -4

Estimated using BJ from data set z

Loss function 2.90859e-021 and FPE 3.18849e-021

Sampling interval: 1 Created: 26-May-2014 15:15:17 Last modified: 26-May-2014 15:15:17 Полученные результаты в среде GUI и в командном окне по параметрам модели совпадают, но по параметрам случайного сигнала различия есть, кроме этого, в результатах расчета в командном окне указывается точность определения коэффициентов в уравнении.

Если на вход объекта подается постоянная составляющая и случайный сигнал от генератора (см. рис. 4), то для дальнейшей работы необходимо выбрать интервал времени (например, t=200...500 с), когда реализация выходного сигнала практически ста-

ционарна (рис. 9, а), а также исключить постоянную составляющую (рис. 9, б). В дальнейшем используется указанная выше методика.

Аналогичные результаты получены и в том случае, когда на вход объекта подается суммарный сигнал f1(t) + u (t) (при этом на вход системы подаются случайный сигнал и сигналы с блоков const или constl, а Gain3=1). Дискретные передаточные функции объекта в указанных случаях отличаются незначительно. Переходные процессы в непрерывном и дискретном объектах и в замкнутой системе по различным каналам практически совпадают (рис. 10).

а б

Рис. 7. Окно задания структуры модели (а) и выбора ее порядка (б)

а б

Рис. 8. Окно сравнения выходов моделей (а) и информации о модели (б)

а б

Рис. 9. Входной и1 и выходной у1 сигналы после выбора временного интервала и исключения постоянной составляющей

И ScoptL In в j 1

а а P^P ASS в

12 10 3 6 4 2 п

И ■ ■>■■■. Г... 1.....Я...1

/.......i

о loo am зоо 4оо 900 Time oilset: 0

И Stopel

а ррр й

Е

О 100 Time offset 0

200 300 400 500

1.5 1

0.5 0

ft.......1.........i.......... ..........

1 jj

О 100 Time offset О

200 300 «0 500

б в Рис. 10. Переходные процессы в непрерывном и дискретном объектах (а), в замкнутой системе

по заданию (б) и возмущению (в)

а

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

C(q) = 1 - 0.1883 (+-0) qA-1 - 0.001941 (+-0.01594) qA-2 - 0.1919 (+-0.0192) qA-3 - 0.4754 (+-0.03425) qM -0.1772 (+ -0.01166) qA-5

D(q) = 1 - 1.68 (+-0) qA-1 + 0.768 (+-0) qA-2 - 0.2472 (+-0) qA-3 + 0.1082 (+-0) qA-4 + 0.05453 (+-0) qA-5

F(q) = 1 - 3.08 (+-0.0003396) qA-1 + 3.525 (+0.0009721) qA-2 - 1.893 (+-0.0009529) qA-3 + 0.5613 (+-

Рис. 11. Переходные процессы в исходной непрерывной и аппроксимированной дискретной системах

по различным каналам

По заданию:

>> z=[y2 u2];th=bj(z,[3 5 5 5 4]); present(th) Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)

B(q) = 0.007084 (+-7.583e-007) qA-4 + 0.02381 (+-1.305e-005) qA-5 - 0.002677 (+-4.215e-005) qA-6

C(q) = 1 - 0.7817 (+-0.1435) qA-1 + 1.219 (+-0.1553) qA-2 - 0.4769 (+-0.2264) qA-3 + 0.6488 (+-0.1322) qA-4 + 0.1571 (+-0.1223) qA-5

D(q) = 1 - 2.127 (+-0.1362) qA-1 + 0.7135 (+-0.36) qA-2 + 0.8842 (+-0.2861) qA-3 - 0.3485 (+-0.03216) qA-4 -0.1212 (+ -0.03048) qA-5

F(q) = 1 - 2.068 (+-0.001764) qA-1 + 1.401 (+0.003492) qA-2 - 0.4107 (+-0.002206) qA-3 + 0.1014 (+0.0006245) qA -4 + 0.004964 (+-0.0001995) qA-5 Estimated using BJ from data set z Loss function 1.60671e-009 and FPE 1.80415e-009 Sampling interval: 1 Created: 27-May-2014 11:44:24 Last modified: 27-May-2014 11:44:24 По возмущению: z=[y2 u2];

>> th=bj(z,[4 5 5 5 4]); present(th) Discrete-time IDPOLY model: y(t) = [B(q)/F(q)]u(t) + [C(q)/D(q)]e(t)

B(q) = 0.002026 (+-2.818e-007) qA-4 + 0.004705 (+-8.824e-007) qA-5

- 0.007729 (+-1.944e-006) qA-6 + 0.0009983 (+-2.514e-006) qA-7

0.000343) qA-4 - 0.113 (+-3.731e-005) qA-5 Estimated using BJ from data set z Loss function 1.72052e-010 and FPE 1.94748e-010 Sampling interval: 1 Created: 27-May-2014 12:59:09 Last modified: 27-May-2014 12:59:10 Исследования на модели показывают, что получение динамических характеристик как объекта, так и системы автоматического регулирования при наличии возмущений в виде случайного сигнала возможно. Однако при распространении метода на промышленные объекты необходимо соблюдать ряд условий.

1. Следует априори установить характер связи между временными рядами данных, так как статистическая связь часто не соответствует причинно-следственной связи.

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

3. В случае использования передаточной функции замкнутой системы для определения динамики объекта необходимо знать применяемый закон регулирования и его параметры настройки.

4. Для определения динамики объекта необходимо получать входную информацию на основании измерения потоков после регулирующего органа.

Статья поступила 03.07.2014 г.

Библиографический список

1. Ротач В.Я. Теория автоматического управления: учебник для вузов. 3-е изд., стереотип. М.: Изд-во МЭИ, 2005. 400 с.

2. Дьяконов В.П. MATLAB 6.5 SP1/7 + Simulink 5/6 в математике и моделировании. Серия «Библиотека профессионала». М.: СОЛОН-Пресс, 2005. 576 с.

3. Баев А.В., Мельник А.А. Применение пакета System Identification Toolbox матричной лаборатории MatLab + Simulink для построения моделей на примере паровых котлов // Вестник Иркутского государственного технического университета. 2013. № 12(83). С. 240-244.