CC BY
78
КОМПЬЮТЕРНОЕ МОДЕЛИРОВАНИЕ
УДК 519.8:666.941
ПОЛУЧЕНИЕ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ОБЖИГА КЛИНКЕРА С ПРИМЕНЕНИЕМ СТАТИСТИЧЕСКИХ МЕТОДОВ
Рассматривается идентификация процессов обжига клинкера с использованием естественных сигналов по различным каналам управления. При этом определяются корреляционные и взаимные корреляционные функции и решается уравнения Винера-Хопфа во временной области. Показана возможность определения структуры сложного объекта управления на основе его декомпозиции на отдельные зоны.
Ключевые слова: идентификация, обжиг клинкера, случайный процесс, корреляционная функция, аппроксимация, переходная функция, комплексная частотная характеристика, система, структура, управление.
Идентификация сложного объекта управления, представляющего собой вращающуюся печь обжига клинкера, с целью дальнейшего использования математической модели для построения системы управления процессом обжига, традиционными способами, основанными на применении типовых сигналов, является невозможной в силу ограничений технологического характера. В связи с этим предлагается использовать для идентификации в качестве исходной информации естественные сигналы по различным каналам управления.
Поскольку существует достаточно большое число управляющих и управляемых переменных, то для иллюстрации подхода к построению динамической модели объекта управления рассмотрим более подробно один из каналов, например, «температура материала в зоне спекания - температура материала в зоне кальцинирования». Статистические данные работы цементной печи для разных каналов были сняты в ОАО «Осколцемент». Как показал анализ случайных процессов, пример реализаций которых приведен на рис. 1, исследуемые процессы относятся к классу эргоди-ческих случайных стационарных процессов в соответствии с определением стационарности в широком смысле.
В.Г. РУБАНОВ В. А. ПОРХАЛО
Белгородский государственный технологический университет им. В. Г. Шухова
e-mail:
Porhalo@mail.ru
Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1
Рис. 1. Реализации случайных процессов: а) температура материала в зоне подогрева, б) температура материала в зоне кальцинирования
Сущность метода получения динамических характеристик объекта состоит в том, что если в нормальных условиях работы воздействие, приложенное ко входу исследуемой системы, может рассматриваться как стационарная случайная функция, то применение коррелятора дает возможность определить корреляционную функцию RQ(т) входа Q(т) и взаимную корреляционную функцию RQq(т) между входом Q(т) и выходом q(т).
Известно [і], что эти функции связаны выражением:
R
Qq (т) = Irq (т - ¿)- w(Ä)dÄ,
(1)
решая которое относительно ю(Х), можем найти импульсную переходную функцию w(t), а по ней и передаточную функцию W(s).
Таким образом, задача определения динамических характеристик объекта сводится к следующим этапам:
- запись случайных процессов на входе и выходе объекта;
- определение вероятностных характеристик случайных сигналов;
- определение динамических характеристик объекта w(t) и W(s). Определение корреляционной функции и взаимной корреляционной функции
входного и въходного сигналов. Исходные формулы для определения корреляционной функции Rq(t) и взаимной корреляционной функции Rqq(t) по экспериментальным реализациям случайных функций имеют вид1:
1
rq (т)=&2T IQ(t)'Q(t+T')dt ’
1 Т
RQq (т) = Tm 2T IQ(t)- q(t + T')dt'
(2)
(3)
Анализируя выражения (2) и (3), можно сказать, что при вычислении RQ(т) и RQq(т) на вычислительной машине, случайные процессы задаются в виде таблиц и в этом случае формулы (2) и (3) необходимо преобразовать.
В результате обработки экспериментальных данных находят не корреляционные функции и мат.ожидания, а их оценки, значения которых вычисляются на ЭВМ по следующим формулам:
0
-Т
Ёа
mr
=0
n + 1
n-k
m
2
i=0
n + 1
RQ (JcAt) = —2(Qi - mQ \Qi+k - mQ
RQq (kAr) =
n - k i=i 1
n - k i=1
2 (Qi- mQ bi
i+k - mq
(4)
(5)
(6)
где Qi=Q[i], qi=q[i] - значения случайного сигнала в i-ом сечении, mQ* и mq* -оценки мат.ожиданий, а Rq* и RQq* - оценки корреляционных функций случайных процессов Q[n] и q[n] соответственно, n - дискретное время. Шаг дискретизации равен 1 мин.
Выражения (4) - (6) определяют алгоритм работы программы на ЭВМ, с помощью которой происходит расчет приведенных характеристик и величин. Программа была написана в среде Borland Delphi 7.0 в виде программного приложения Windows. Исходные данные находятся в текстовых файлах. Также была проведена проверка с помощью аналогичных вычислений в среде Microsoft Excel 2003. Полученные результаты совпали.
В результате вычислений были получены значения для mQ* и mq*, равные соответственно 137.7 0С и 245.5 0C, а также значения оценок корреляционной функции Rq*(t) и взаимной корреляционной функции RQq*(r), которые были нормированы. Графики функций Rq*(t) и RQq*(T) и их аппроксимация приведены на рис.2 и рис.3 соответственно.
Rc
1.2
0.6
0.2
0 0 5 )
И
Рис. 2. Оценка корреляционной функции для температуры материала в зоне подогрева
::
Как видно из приведенного графика, вычисленная оценка мат. ожидания т^‘ является достаточно точной, поскольку с течением времени оценка корреляционной функции стремится к нулю.
Следует учитывать, что расчетные значения имеют не вполне плавный характер изменения, однако тенденция поведения корреляционных функций угадывается отчетливо, что позволяет аппроксимировать их плавными кривыми. Вычисленные оценки корреляционных функций отражают характер рассмотренных стационарных процессов. Действительно, относительно медленное затухание характеристик (20 минут) указывает на то, что они должны соответствовать случайным процессам, изменяющимся медленно, что вполне соблюдается.
Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1
Рис.3. Оценка взаимной корреляционной функции между температурой материала в зоне кальцинирования и температурой материала в зоне подогрева
Определение импульсной переходной функции Шо(Ь) и переходной функции Решение интегрального уравнения Винера-Хопфа (1) аналитическими методами требует представления функций и RQq* в виде аналитических выражений. Сложность аналитических выражений зачастую приводит к неразрешимости данного уравнения в силу нахождения искомого в подынтегральном выражении, а переход в частотную область требует применения преобразования Фурье. Так как в расчетах значения оценок корреляционных функций заданы в виде последовательности ординат, следующих друг за другом через дискретные отрезки времени Ат, то решение интегрального уравнения (1) удобно свести к решению системы алгебраических уравнений с числом уравнений и неизвестных, равным количеству вычисленных ординат статистических характеристик. Для этого интеграл в уравнении Винера-Хопфа аппроксимируется конечной суммой:
N
1= У w\
(кАт) = ^ w(kAт) ■ Я*д (кАт - пАт)Ат,
п= 0
разделив выражение (7) на Ат, получим:
N
Яп ^(пАт)■ Квп ’
(7)
(8)
п= 0
Я0п (¡Ат) *, ч
где qi =—-, Юи=ю(п Ат ) и Я0п = Я0(кАт - пАт).
Ат
Равенство (8) можно представить в виде системы п линейных уравнений с п неизвестными:
'^11^ + ^012^ + К0^3 + К + &0^п = ql >
К0 21 wl + К022 ^2 + ^ 23 W3 + К К0 2nwn = q2’ )
Квп№ + КОп2 W2 + КОп3 W3 + к + Явпп^ = Яп ■
В случае вычисления корреляционных функций нами было получено 21 значение для ординат, поэтому для нахождения ординат функции ю(і), необходимо ре-
шить систему уравнений (9) при п=21. Применение ЭВМ для решения данной системы не представляет особых затруднений.
Однако следует помнить, что система уравнений (9) является плохо обусловленной и даже малые изменения значений ординат оценок корреляционных функций могут приводить к значительным скачкам искомых значений ординат весовой характеристики ю(Ь). Решение системы (9) можно дать, воспользовавшись представлением в матричной форме:
Ш = Л~^, (10)
где Л-1 матрица обратная А. Матрица существует только для квадратной несобственной матрицы, т.е. если определитель прямой квадратной матрицы А не равен нулю.
Выражения для отдельных матриц имеют вид:
А =
Яв (0) Яв (Ат) Яв (2 Ат)
Кв (Ат) ........ Яв [(п - 1)Ат]
Яв (0)
Яв [(п - 2)Ат]
Яв (Ат) ....... Яв [(п - 3)Ат]
Яв[(п - 1)Ат] Яв[(п - 2)Ат]
Яв (0)
(11)
Ж = ш1 w1
w„
в = ІІЯі Яі
Яп
(12)
Решение системы (9) производилось в пакете MathCAD 2001. Кривая, построенная по полученным результатам, приведена на рисунке 4,а.
Рис. 4. Импульсная переходная и переходная функции
Известно [2], что переходная функция представляет собой интеграл от импульсной переходной функции к^(Ь).
Для определения переходного процесса в системе вызванного произвольным воздействием /(і), в случае, когда известна весовая характеристика, можно воспользоваться уравнением свертки:
х'в(і) = |wQ (і-т)/(т)іт .
(13)
Переходная функция это реакция на единичное воздействие, тогда, учитывая, что/(Ь)=1 и, заменяя Ь-т = X, получим:
кв (?) = |кв (А)йй .
(14)
0
Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1
Так как в результате численного метода решения уравнения (9) весовая характеристика получена в форме графика, то для нахождения переходной функции был применен пакета МаШСАБ 2001(рис.4,б).
Анализ переходной характеристики показывает, что при сравнительно быстром изменении температуры материала в зоне подогрева, температура материала в зоне кальцинирования сначала растет, а затем происходит медленное ее выравнивание относительно нового состояния, что отвечает реальным условиям течения процесса.
Расчет комплексной частотной характеристики объекта по кривой переходного процесса. Если известна переходная характеристика И.р@), заданную графически, то для нахождения комплексной частотной характеристики (КЧХ) необходимо заменить кривую переходного процесса ломаной линией, т.е. аппроксимировать ее кривой Нд (і), состоящей из сопрягающихся друг с другом прямолинейных отрезков.
Тогда Нд (і) = -д (і) .
К(і)
і
п
Г
і
Є.
I, МНИ ->
Рис. 5. Аппроксимация переходной функции hQ(t) ломаной линией
Построенная кривая Нд (і), отражает переходную функцию hQ(t), как показано
на рис.5. Определив, далее коэффициенты, характеризующие углы наклона полученных прямолинейных отрезков:
Нд1 Нд г Нд1
к— = -д—; к2 = -д2--^; ... (15)
¿і Ь - і—
можно найти значение вещественной и(ш) и мнимой У(ш) частей комплексной частотной характеристики W(jы)[з]:
і
и(ю) = — У кі (біп ґію - БШ іі-ію),
0 і=і
1 п
V(ю) = — У кі (собі{ю- СОБіі-ію).
і=і
Введя обозначение: ті = — — 1 ; Аі = ———; <жі = 0і; 2кіті = Ві; — = уі, и
2 2 Ті
осуществив тригонометрические преобразования, получим расчетные формулы для определения частотных характеристик:
и (<) = Х Б ^ Оі),
і =1
V (<) = -£ Б £ (V, О,).
(17)
=1
В соответствие с изложенной выше методикой получим график КЧХ для объекта, имеющего переходную характеристику
График КЧХ Ш(]ш) представлен на рис. 6. По виду комплексной частотной характеристики объекта исследования можно определить передаточную функцию, как динамическую модель объекта в области комплексного переменного.
т<ч
0.4 ■ 0.5 о) = 1 )
со—ЗС со-0
2 Д 0 2 0 4 0 6 0 со 8 =0.05 1 2 1 4 6 і
і
\ (0=3
V )=10
\
со—5
1.6 ■
I :<(■>)
Рис. 6. Комплексная частотная характеристика объекта
Модель печи в целом может быть представлена в виде распределенных в пространстве, примыкающих друг к другу зон, преобразования материала. В этом случае целесообразно разделить печь на п зон по количеству контролируемых переменных. На каждом участке последовательной цепи моделей происходят физико-химические превращения материала, отражаемые модельными переменными в форме оператора преобразования, причем контролируемая переменная на выходе предыдущей зоны влияет на последующие зоны печи вплоть до выхода объекта. Такой подход к математическому описанию сложного объекта управления существенно упрощает теоретическое описание исследуемого объекта, поскольку декомпозиция его на более простые объекты позволяет глубже проанализировать реальные процессы в условиях, присущих конкретной зоне.
Располагая реальными случайными процессами, описывающими характер изменения температуры по зонам печи, изменением подачи топлива во времени, а также статистическими характеристиками, полученными в соответствии с изложенной теорией, устанавливающими связь между температурой в зоне и содержанием окиси кальция в смеси, можно представить обжиговую печь в форме структурной модели (рис. 7), в которой структурными блоками являются передаточные функции каждой из зон печи по каналу «температура материала - содержание свободной окиси
Серия История. Политология. Экономика. Информатика. 2010. № 7 (78). Выпуск 14/1
кальция» - УгСв) и каналу «количество теплоты на обжиг - температура материала» -Wi(s), а в качестве переменных выступают физико-химические свойства материала Х(Ь), которые проявляются в г-ой зоне и воспринимаются в форме контролируемой переменной формируемой прибором первичной информации (температура материала); рг(Ь) - количество тепла, подводимое в г-ю зону; аг(Ь) - входная переменная г-ой зоны, характеризующая физико-химические свойства материала, например содержание свободной окиси кальция на выходе зоны; д$) - аддитивная помеха, х$) -переменная используемая для синтеза управления; у@) - выходная переменная печи.
Рис. 7. Структурная схема модели процесса обжига Литература
1. Рубанов, В. Г. Статистическая динамика систем управления [Текст] : учеб. пособие / В.Г. Рубанов. - Белгород: БелГТАСМ, 2000. - 113 с.
2. Бессекерский В. А., Попов Е.П. Теория систем автоматического управления [Текст] / В.А. Бессекерский, Е. П. Попов - Спб: Профессия, 2003. - 757 с.
3. Фатеев А. В. Основы линейной теории автоматического регулирования [Текст] / А.В. Фатеев - Л.: Госэнергоиздат, 1954г. - 296 с.
THE MATHEMATICAL MODEL IDENTIFICATION OF A BURNING OF THE CEMENT CLINKER, WITH THE STATISTICAL METHODS USAGE
V. G. RUBANOV V. A. PORKHALO
Belgorod State Technological University named after V.G. Shoukhov
There considered the identification of the process clinker burning using source signals from different control channels. For this, calculated correlation and mutual correlation functions and the Wiener-Hopf equations solved in time plan.
Key words: identification, clinker, random process, correlation function, approximation, the transition function, the complex frequency response, system, structure, control.
e-mail: Porhalo@mail.ru
