Structural and parametric identification of the process model using a rotary pulsation machine.pdf Проблема идентификации моделей динамических объектов и технологических процессов является одной из актуальных в современной мировой науке и практике [1-16]. Математические модели используются практически во всех отраслях современной промышленности, и именно с их применением связаны успехи систем автоматизации действующих производств. Не является исключением и пивоваренная индустрия, которая за последние десятилетия осуществила прирост производственных мощностей. Начало ХХІ в. ознаменовано стремительным техническим развитием отечественной пивоваренной отрасли, направленным на увеличение объемов производства пива, расширение его ассортимента и повышение качества. В пивоварении технологический процесс охмеления пивного сусла играет важную роль. В процессе охмеления происходит ряд технологических явлений, таких как экстрагирование и изомеризация горьких веществ хмеля, придание хмелевой горечи, ароматизация пивного сусла и т.д. Применяемые в промышленности способы охмеления пивного сусла не дают достаточно полного и эффективного выхода экстрактивных веществ, специфических компонентов хмеля. Поэтому в последнее время в российской промышленности одним из наиболее перспективных способов интенсификации процесса охмеления пивного сусла является использование роторно-пульсационных аппаратов (РПА) [17-22], которые позволяют максимально извлекать горькие вещества из хмеля. Для изучения закономерностей влияния технических и технологических характеристик РПА на процесс охмеления пивного сусла представляется актуальным проведение математического моделирования этого технологического процесса. Первые попытки моделирования процессов на выходе РПА были осуществлены на основе регрессионных моделей планирования эксперимента [23], однако для их построения необходим большой объем экспериментальных данных, что не всегда является практически реализуемым. Также имеется ряд исследований по определению влияния гидромеханических и акустических явлений в РПА на характеристики получаемого продукта [17-19]. В работах [20-22] проводится численное моделирование процессов растворения в данных аппаратах. Однако на сегодняшний день вопрос моделирования динамических процессов в РПА остается открытым. В связи с этим представляется актуальным использование альтернативных методов моделирования данного технологического процесса. Одним из таковых является кибернетический подход к анализу РПА. В работе на основе алгоритмического аппарата непрерывных дробей будет рассмотрена структурно-параметрическая идентификация, позволяющая построить динамическую модель получения хмелевого экстракта в условиях отсутствия априорной информации и при минимальном объеме экспериментальных данных. 63 М.А. Новосельцева, С.Г. Гутова, Е.С. Каган, Д.М. Бородулин 1. Постановка задачи структурно-параметрической идентификации динамического объекта Априорная неопределенность и меняющиеся условия функционирования динамических объектов являются характерной чертой научных исследований, что значительно затрудняет применение большого количества существующих методов идентификации [1-17], в которых восстановление структуры модели является неочевидным процессом и приводит к перебору пробных моделей на основе субъективного экспертного подхода. Поэтому с целью принятия научно обоснованных технических решений становится актуальным использование методов структурно-параметрической идентификации. Линейный динамический объект описывается непрерывной передаточной функцией (НПФ) Y (J) b,S' + b.s \\ / m_m-1 m-1 -b2 - bs + bg a„s ~^,s" 1 +... + as2 + as + (1) где s - переменная преобразования Лапласа [3]. Известно также модальное представление НПФ объекта [3, 24]: G(s) = K n(s - s0 yn(s - s,), (2) i=1 / i=1 где K = bmlon, si0, S20, ^, Sm0 - нули НПФ, Si, S2, ^, Sn - полюсы НПФ. Использование современных систем управления и новейших информационных технологий приводит к цифровым процессам получения и обработки информации, поэтому входное x(t) и выходное y(t) воздействия объекта могут быть представлены результатами измерений с некоторым временным шагом А: x[0] = x(0A), x[1] = x(1A), x[2] = x(2A), _ y[0] = y(0A), y[1] = y(1A), y[2] = y(2A)... Тогда модель динамического объекта в цифровой форме с использованием преобразования z = esA [24] имеет вид: (3) bm^ + ... + b1z b az" +...+az1 +1 где G(z) - дискретная передаточная функция (ДПФ) динамического объекта [24, 25], X(z), Y(z) -z-преобразования числовых последовательностей значений входного и выходного воздействий. Дискретную модель объекта можно представить и в форме конечно-разностного уравнения у (nAt )=Z b,x[n - i ]-2 [n - i]. (4) i=0 i=1 Таким образом, задачей структурно-параметрической идентификации динамического объекта является построение моделей вида (1)-(4) по результатам измерений входного и выходного воздействий с шагом дискретизации A. 2. Структурно-параметрическая идентификация динамического объекта на основе теории непрерывных дробей В последние десятилетия аппарат непрерывных дробей достаточно активно развивается в области фундаментальных исследований их свойств [26-28]. Однако вопросу использования непрерывных дробей в прикладных исследованиях отводится очень малая роль. Ранее в авторских работах [24, 25, 29] теория непрерывных дробей была успешно использована для решения широкого спектра прикладных задач. Таковыми являлись задачи структурно -параметрической идентификации моделей различных динамических объектов, анализа свойств и характеристик детерминированных, периодических и стохастических процессов, оценки риска объектов и систем и т.д. Для построения модели объекта поместим результаты измерений входного и выходного воздействий в строки № 0 и № 1 матрицы-идентификатора (5): 64 Структурно-параметрическая идентификация модели технологического процесса столбец № 0 столбец №1 столбец № 2 столбец строка № 0 X [0] X[1] X[2] X[3] строка № 1 У [0] У [1] У [2] У [3] строка № 2 a20 a21 a22 a23 строка № 3 a30 a3i a32 a33 (5) (6) Остальные элементы матрицы-идентификатора находятся по формулам am,k ~ am-2,k+1 am-2, 0 “ am-\\,k+\\lam-\\, 0 , ao,k = x[k], ai,k = y[k], k = 0,1,2,... Расчет матрицы-идентификатора прекращается при появлении строки с нулевыми элементами. Данный факт является критерием останова вычислительной процедуры по формуле (6). Элементы столбца № 0 матрицы (5) образуют непрерывную правильную С-дробь [24, 25], аппроксимирующую ДПФ динамического объекта: aio Zaoo G(z ) = ■ 1 + a20z 1 (7) В случае j[0] = 0 строка № 1 матрицы сдвигается на один элемент влево. Отрицательная степень переменной z для данной строки увеличивается на единицу. Затем непрерывную дробь (7) сворачивают в дробно-рациональное выражение, аналогичное (3), и определяют ее нули и полюсы. Для восстановления модели объекта в форме НПФ (2) требуется использовать соотношение 5 = (іП z| + iargz^A. (8) При усложнении объекта у ДПФ могут появляться дополнительные нули и полюсы [24, 25], которые отсутствуют в НПФ. Это отрицательные или равные нулю корни числителя и знаменателя (7). Дополнительные нули и полюсы по формуле (8) в 5-плоскость не переводятся [24, 25]. Далее находят коэффициент K согласно теореме о конечном значении [3] основных свойств z-преобразования: lim G ( z ) = K. (9) Пример. В качестве иллюстрации работы метода структурно-параметрической идентификации рассмотрим апериодический объект 1-го порядка. Его НПФ имеет вид [3]: G(i ) = (10) K Is + 1 где T - постоянная времени, K - коэффициент усиления. Реакция объекта на единичное ступенчатое воздействие описывается разгонной характеристики h(t) вида: h(t ) = K (і - е). (11) Зададим некоторый шаг дискретизации А и заполним матрицу-идентификатор: 1 1 1 K (1 - еА 'т) K (1 - e-^A/I) K (1 - e-3 a20 a21 a22 0 0 0 (12) Строка № 0 определяется значениями входного единичного ступенчатого воздействия. Элементы строки № 1 - значения разгонной характеристики (11), рассчитанные с шагом дискретизации А. Элементы строки № 1 сдвинуты на один столбец влево. Элементы строки № 2 рассчитаны по формуле (6): 65 М.А. Новосельцева, С.Г. Гутова, Е.С. Каган, Д.М. Бородулин /«00 - «u.i /«10 = ■ (1 - ) / (1 - ) причем -д/г Элементы строки № 3 равны азк = 0. Появление строки с нулевыми значениями определяет ДПФ объекта: a0ila00 ■_ K(1 - е^д/г )z ^01l ^-*00 1 ^ a20 z G (z) = : В ДПФ (13) нулей нет, полюс равен что при переходе в 5-плоскость дает 51 = -HT. Последний результат совпадает с полюсом объекта (10). Тогда модель динамического объекта (10) имеет вид: kT 1 -д/г -1 1 - e z zi = G(s ) = k-1- (13) s +1 /г Ts +1 Для определения коэффициента передачи kT воспользуемся формулой (9) и получим lim G ( z) = K . z ^1 Дискретная модель разгонной характеристики принимает вид: h[n] = е~д/т • h[n-1]+K(1 -е~д/T)■ x[n-1]. (14) Рис. 1. Разгонная характеристика (11) и значения ее дискретной модели (14) (сплошная черная линия -график непрерывной разгонной характеристики, синие квадрата: - значения дискретной модели) Fig. 1. Overclocking characteristic (11) and values of its discrete model (14) (solid black line -graph of continuous acceleration characteristics, blue squares-values of discrete model) На рис. 1 для сравнительного анализа приведены графики значений разгонной характеристики (11) и ее дискретной модели (14) для А = 1, K = 2, T = 1. Очевидно, что дискретная модель точно оценивает значения истинной разгонной характеристики объекта и может быть использована для оценки состояния процесса на выходе объекта. Таким образом, на основе цифровых отсчетов входного и выходного воздействий апериодического звена 1 -го порядка с использованием метода структурнопараметрической идентификации были определены структура и параметры модели объекта (10). 3. Структурно-параметрическая идентификация модели технологического процесса получения хмелевого экстракта Разработка математической модели, описывающей технологический процесс получения хмелевого экстракта на выходе РПА, является актуальной задачей для пивоваренных отраслей промыш-66 Структурно-параметрическая идентификация модели технологического процесса ленности, поскольку дает возможность определить необходимые технологические параметры работы (температура, частота вращения ротора), которые позволят максимально улучшить показатели по выходу экстракта. Осуществим структурно-параметрическую идентификацию этого процесса с использованием описанного выше метода. Входным воздействием на РПА является количественное соотношение неохмеленного сусла и гранулированного хмеля F. Воздействие на выходе РПА - показатель эффективности процесса - содержание изогумулона в хмелевом экстракте (в мг/л). Техническими и технологическими параметрами, влияющими на процесс охмеления пивного сусла хмелевыми экстрактами, являются температура перерабатываемой среды t (°С), частота вращения ротора n (об/мин) и величина зазора между ротором и статором, измеряемая в миллиметрах. С помощью лабораторных исследований экспертами было установлено, что при величине зазора между ротором и статором, равной 0,3 мм, наблюдается лучший выход изогумулона в хмелевом экстракте. Поэтому при проведении дальнейших расчетов будет использоваться именно это значение зазора. Зафиксируем значения параметров t = 55°С, n = 2 500 об/мин. Экспертами были предоставлены значения разгонной характеристики1 по содержанию изогумулона, измеренные с шагом А = 1 мин, в качестве входного воздействия было взято постоянное соотношение неохмеленного сусла и гранулированного хмеля, равное 8 мг/л. Рассчитаем матрицу-идентификатор, в которой строка № 0 занята значениями входного воздействия, а строка № 1 - значениями разгонной характеристики: 8 22, 40 8 18,82 -0,065887 -1,764758 1,876144 8 20,06 -0,186504 -1,684464 1,908832 8 22,33 -0,189161 -1,190223 8 22,38 -0,190223 В строке № 1 осуществлен сдвиг на один столбец влево, так как значение разгонной характеристики в начальный момент времени равно 0. Значения последующих строк рассчитаны по формуле (6). Строка матрицы № 5 может считаться равной нулю с точностью проведенных измерений. Элементы столбца № 0 порождают непрерывную дробь вида: G(z ) = ■ 1 - 18,82 -1 -Z 8 0,065887z -1 1,764758z -1 + 1,876144z Сворачивая полученную непрерывную дробь, получим дробно-рациональное выражение для ДПФ объекта . S 2,3525z-^ + 0,262035z ^^) ” 1 + 0,045499z-' - 0,123614z ^ 1- (15) В 5-плоскость переводится только один полюс - zi = 0,329574, получаем si = -1,109. Тогда T = 0,9009, K = 2,836078. НПФ примет вид: 2,8361 G(j) = 0.9009s +1 (16) а конечно-разностное уравнение h [n] = -0,045499 • h [n -1] + 0,123614 • h [n - 2] + 2,3525 • x [n -1] + 0,262035 • x [n - 2]. (17) 1 Данн^іе собран^! кандидатом технических наук, доцентом кафедры «Технологическое проектирование пищевых произ-водств2 Института инженерн^іх технологий Кемеровского государственного университета Е.А. Сафоновой. 67 М.А. Новосельцева, С.Г. Гутова, Е.С. Каган, Д.М. Бородулин На рис. 2 приведена линейная аппроксимация исходных данных и полученной дискретной модели. Результаты идентификации при других значениях температур и частотах вращения ротора приведены в табл. 1, 2. Рис. 2. Линейная аппроксимация экспериментальн^іх данных и полученной дискретной модели разгонной характеристики (сплошная черная линия - значения конечно-разностного уравнения, серая линия с квадратами - линейная аппроксимация экспериментальн^іх данн^іх) Fig. 2. Linear approximation of the experimental data and the resulting discrete model of the acceleration characteristic (solid black line - values of the finite difference equation, gray line with squares-linear approximation of the experimental data) Таблица 1 Значения K при различных t и n П1 = 2,5 П2 = 3 Температура t, °С t1 = 55 2,836078 3,388565 t2 = 70 3,209291 3,879953 t3 = 85 3,605610 4,279379 Значения T при различных t и n T Частота вращения ротора n (1000 об/мин) П1 = 2,5 П2 = 3 Температура t, °С ?1 = 55 0,9009 0,9950 t2 = 70 0,7348 0,8292 t3 = 85 0,7610 0,7407 Таблица 2 Сравним результаты идентификации процесса охмеления с результатами идентификации апериодического звена 1-го порядка. Наличие дополнительных нуля и полюса ДПФ в модели процесса охмеления позволяет сделать предположение о том, что объект идентификации имеет более сложную структуру и не является линейным. На основе полученных результатов зададим зависимости K = K(t, n) и Т = Т (t, n) и далее рассчитаем их значения с помощью метода наименьших квадратов: K = -0,6886 + 0,0126? + 0,6986и + 0,0081/и - 0,0001t2 , T = 0,9340 - 0,0276? + 0,6461n - 0,0076?n + 0,0003?2. (18) Ошибки моделирования, найденные с помощью квадратического критерия, оказались равны: ^[Kij - k( t, Uj)) 2= 0,001, ^{Tij - t( t, nj)) 2= 0,004. Таким образом, для значений параметров t, n (55 < t < 85, 2,5 < n < 3) модель процесса охмеления пивного сусла на выходе РПА описывается НПФ вида: K (?, п) T (?, п) s +1 ’ 68 Структурно-параметрическая идентификация модели технологического процесса а разгонная характеристика - следующей формулой: (19) h(t )=K

Скачать электронную версию публикации