CC BY
98
УДК 66.012-52
МОДЕЛИРОВАНИЕ ПРОМЫШЛЕННЫХ ОБЪЕКТОВ УПРАВЛЕНИЯ С ДЕТЕРМИНИРОВАННОЙ НЕСТАНЦИОНАРНОСТЬЮ
В.Ф. Беккер, Н.В. Бильфельд, Д.С. Камаев*
Пермский национальный исследовательский политехнический университет, Березниковский филиал *Филиал «Ависма» ОАО «Корпорация ВСМПО-Ависма», г. Березники E-mail: [email protected]
Рассматривается актуальная задача исследования динамики нестационарных промышленных объектов средствами и методами математического моделирования. При использовании стандартных средств моделирования из библиотеки Simulink отсутствует возможность изменения параметров моделируемого звена во времени в процессе решения задачи. Для устранения этого недостатка предложен способ моделирования промышленных объектов с динамически меняющимися параметрами на моделях, содержащих автономный модуль, реализующий произвольную функциональную зависимость параметров моделируемых стандартных звеньев от машинного масштаба времени. В основу моделирования положена реализация интегрирующего звена с параметром, изменяющимся во времени по определенному закону, или в зависимости от внешнего сигнала любой формы. Показан переход к моделированию инерционного звена первого, второго и более высоких порядков. Приведен пример двухпозиционной релейной стабилизации температуры воды в баке при различных внешних возмущениях в проточном реакторе идеального перемешивания. Обсуждаются результаты вычислительного эксперимента. Выявлены ограничения на области применения модели.
Ключевые слова:
Детерминированная нестационарность, математическое моделирование, двухпозиционное регулирование, среда Simulink.
Промышленные объекты управления часто являются существенно нестационарными. Влияние нестационарности объекта управления на динамику системы иногда удается описать детерминированной зависимостью, которая на основе фундаментальных законов переноса массы и энергии позволяет определить скорость протекания технологического процесса, а следовательно, и кинетическое время, необходимое для достижения конечного состояния системы или завершенности процесса при данной скорости. Примером такой нестацио-нарности является промышленный аппарат вакуумной дистилляции губчатого титана, получаемого магнийтермическим восстановлением из тетрахлорида (Кролл-процесс) [1]. В этом аппарате скорость высокотемпературной отгонки примесей пропорциональна изменяющемуся во времени количеству примесей, находящихся внутри аппарата. Другим примером такой детерминированной нестационарности является высокотемпературная отгонка под вакуумом паров цинка [2].
Рассмотрим задачу исследования динамики нестационарных промышленных объектов, которые практически всегда с приемлемой для практики точностью удается описать последовательным включением нескольких инерционных звеньев первого порядка с запаздыванием [3, 4]. Выполняя последовательную декомпозицию такого объекта, получаем элементарное звено вида W(s)=k/(Ts+1).
При моделировании в среде Simulink параметры звеньев задаются в меню Parameters. Запуская затем процесс решения модели, нужно учитывать, что в ходе вычислительного эксперимента эти параметры (в частности постоянная времени T) сохраняют введенные изначально значения, так как в стандартном наборе средств Simulink отсутствует возможность изменения параметров моделируемого звена во времени. Как при моделирова-
нии в таких жестких условиях изменять свойства моделируемого объекта в процессе проведения вычислительного эксперимента?
Выход был найден с использованием стандартных блоков библиотеки БтиШк. Рассмотрим реализацию интегрирующего звена вида 1/Тэ с переменным параметром Т в виде произведения
— =------, в котором множитель 1/Т можно менять
Ts Т s
динамически, в процессе работы модели (рис. 1).
Рис. 1. Реализация в Simulink. интегратора с переменным параметром T
Здесь основной сигнал поступает с блока Step и представляет собой единичный скачок в момент времени t=1 с. Блок Pulse Generator задаёт изменение параметра T с периодом 2 с и амплитудой 4, а сложение с константой 1 поставлено для того, чтобы избежать деления на ноль. То есть параметр T меняется от 1 до 5. Таким образом, нарастание сигнала на выходе интегратора будет обратно про-
порционально величине динамически изменяющегося во времени параметра Т. Полученные графики представлены на рис. 2.
Важным моментом является то, что параметр Т интегратора можно менять не только по заранее определённому закону, но и, например, в зависимости от выходного сигнала интегратора или любого другого блока реализуемой модели. Единственное условие, как уже отмечалось выше, нужно предусмотреть (исключить) возможность деления на ноль.
Рис. 2. Зависимость выходного сигнала интегратора от параметра Т
Теперь рассмотрим интегрирующее звено, но с отрицательной обратной связью с коэффициентом к (рис. 3). Известно, что такой интегратор эквивалентен инерционному звену первого порядка. Покажем это.
Ж (5)
Wos (s) =
1 + W (sWos (s) 1
Wos (s) =
1 + -k - (s + к ) к i - s +1 s s \ к )
= TW-(s),
= T-
1
Ts +1
где Т=1/к, а ^(я) - передаточная функция инерционного звена первого порядка.
Таким образом, чтобы выразить ^(з) через WOS(s), достаточно разделить её на параметр Т. А так как в сочетании с простым интегратором мы можем динамически менять параметр Т, то значит, и для инерционного звена первого порядка можем сделать то же самое.
Рис. 3. Интегратор с отрицательной обратной связью
Удобно организовать такое звено в виде подпрограммы (подсистемы) - блок Subsystem из раздела Ports & Subsystems библиотеки Simulink. Пример такой реализации изображён на рис. 4.
Такая подсистема полностью идентична звену с передаточной функцией
W-( s) = -±-.
Ts +1
Рис. 4. Инерционное звено первого порядка с переменным параметром Т, реализованное в виде модуля
Рассмотрим проточный реактор идеального перемешивания. Как известно из классического моделирования, такой реактор представляется инерционным звеном первого порядка с постоянной времени Т, которая характеризует время пребывания жидкости в реакторе. Однако это справедливо, лишь когда это время действительно постоянно. А если, например, входящий поток имеет значительные колебания, или уровень жидкости (то есть фактически объём реактора) меняется? В этом случае соответственно будет меняться и «постоянная» (а фактически «переменная») времени T реактора. Поэтому будем называть её просто «параметр T» реактора.
Итак, рассмотрим случай, когда входной поток в реактор переменный. Например, бак, служащий для стабилизации температуры воды (рис. 5).
Вода
Вода
20°
w Is)
р /0/
ч /
Рис. 5. Система стабилизации температуры воды в баке
Холодная вода подаётся постоянно, а подачу горячей воды регулирует двухпозиционный регулятор. Цель регулирования - поддержание температуры на выходе бака в пределах 40 ... 60 °С. Объём бака примем 7=200 л. Поток холодной воды 0соо1=2 л/с, а горячей - 0Ь(Л=10 л/с. Исходящий поток будет равен сумме входящих. Температуры потоков примем Тсоо1=20, ТЬо1=80 °С. Модель, реализованная в БтиНпН, представлена на рис. 6.
Рис. 6. Модель системы стабилизации температуры воды в баке
Температура воды после смешивания на входе в бак, принимая постоянной теплоемкость воды, будет равна
T • Q + T • Q
T = cool -g^cool______hot -g^hot
mix _ Q + Q '
cool hot
Параметр T бака будет определяться отношением ^/(Qcooi+Qhot).
Регулятор реализован с помощью реле с гистерезисом. На его выходе формируется «1» для открытия клапана горячего потока при снижении темпе-
ратуры в баке до 40 °С, и «0» для закрытия при достижении 60 °С. Двухпозиционный регулятор здесь взят для наглядности сравнения переходных процессов при разных величинах входящего потока. Сами переходные процессы представлены на рис. 7.
Бак в этой модели представлен подсистемой ІпетИ, в которой реализована схема, изображённая на рис. 4. Начальная температура в баке принята 20 °С.
Из графиков видно, что когда клапан подачи горячей воды открыт, температура в баке изменяется гораздо быстрее, чем при закрытом клапане. Причём, даже не смотря на то, что при открытом клапане начальное отклонение температуры на входе в бак меньше (|70—40|°С < |20-60|°С). Это соответствует уменьшению параметра Т бака при возрастании потока.
Рис 7. Изменение температуры в баке при разных входных потоках
Но применение блоков с переменным параметром не ограничивается инерционными звеньями первого порядка. Рассмотрим реализацию инерционного звена второго порядка с помощью того же подхода.
Звено второго порядка можно представить в виде последовательного соединения двух звеньев первого порядка:
W2(s) =
1
T- s + T2 s +1
(T-s + 1)(T2s +1)
или
W2(s) =
1
1
T s + T2 s +1 TTs + (T + T2)s +1
(i)
(2)
где Ж2(я) - передаточная функция инерционного звена 2-го порядка; Т1* и Т2* - постоянные времени звена второго порядка; Т1 и Т2 - постоянные времени последовательных звеньев первого порядка.
Теперь, в зависимости от того, в каком виде дано исходное звено второго порядка - (1) или (2), можно смоделировать его в виде подсистемы с переменными параметрами.
Пусть звено второго порядка задано в виде (1). Тогда модель с переменными параметрами реализуется простым последовательным соединением. Пример такой подсистемы представлен на рис. 8. Особенностью этого варианта является то, что поскольку при последовательном соединении линейных звеньев их порядок не важен, параметры Т1 и Т2 можно менять местами.
Если же звено второго порядка задано в виде (2), и параметры Т1 и Т2 неизвестны, то они вычисляются из системы уравнений
\Т • Т2 = Т•;
[Т + Т2 = Т2\
которая сводится к квадратному уравнению. Решение квадратного уравнения также можно реализовать в рамках подсистемы. Пример такой реализации представлен на рис. 9.
Особенность данной реализации заключается в том, что при вычислении квадратного корня из дискриминанта блок sqrt не даёт ошибку при отрица-
тельном дискриминанте. То есть если D<0, то sqrt(D)=-sqrt|D|. Поэтому вычисленные параметры Ті и Т2 никогда не будут содержать мнимой части, но при D<0 они будут найдены неверно. Это нужно учитывать при выборе диапазона значений Т1* и Т2*, подаваемых на вход блока.
Рис. 8. Модуль реализации инерционного звена 2-го порядка при известном разложении на последовательное соединение звеньев 1-го порядка
Аналогично можно смоделировать и систему третьего порядка. Поскольку при неизвестном разложении её передаточной функции на множители решение кубического уравнения через дискриминант весьма громоздко, гораздо эффективнее использовать функцию roots из пакета MatLab. Для этого используем блок MATLAB Function из раздела User-Defined Functions библиотеки Simulink. На рис. 10 представлен пример реализации модели объекта с передаточной функцией вида «инер-
Рис. 9. Модуль реализации инерционного звена 2-го порядка при неизвестном разложении на последовательное соединение звеньев 1-го порядка
Рис. 10. Инерционное звено третьего порядка с переменными коэффициентами
ционное звено третьего порядка», когда полином третьей степени записывается в следующем виде:
Щ( 5) = -----------*12----;----- =
3 Т1*53 + Т2*52 + Т3* 5 + 1
=_____________________________1____________________________=
= Т1Т2Т353 + (Т1Т2 + Т1Т3 + Т2Т3)52 + (Т + Т2 + Т3)5 +1 = 1
мнимые. И в дальнейшем переменные параметры подавать только на чисто инерционную часть.
Block Parameters: roots
MATLAB Fen
(Ts + 1)(T2s+ 1)(T3s +1)
Подаём на соответствующие входы модели переменные параметры T1*, T2* и T3*, а также константу полинома (свободный член) на вход «1». В блоке MATLAB Function параметры устанавливаем, как показано на рис. 11. Блок Complex to Real-Imag служит для отделения вещественной части корней полинома от мнимой. Дело в том, что блок Integrator не может работать с комплексными числами. До тех пор, пока все корни полинома вещественные, такая модель работает. Но как только появляется хоть одна пара комплексных корней, процесс моделирования принудительно останавливается, т. к. в них вещественная часть становится неравна самому корню, что ведёт к неверным результатам вычислений.
Пара мнимых корней обычно свидетельствует о наличии колебательной составляющей в объекте. На практике такие объекты встречаются не часто. В этом случае можно рекомендовать попытаться разбить передаточную функцию на две части: в одной только комплексные корни, а в другой только
Pass the input values to a MATLAB function for evaluation. The function must return a single value having the dimensions specified by 'Output dimensions' and 'Collapse 2-D results to 1-D'.
Examples: sin, sin(u], foo(u(1), u(2))
Parameters— MATLAB function:
|ИЩ
Output dimensions:
3
Output signal type: complex Collapse 2-D results to 1 -D
E
OK
Cancel
Help
Apply
Рис. 11. Параметры блока MATLAB Function для звена третьего порядка
В общем случае любой полином, не содержащий комплексных корней, раскладывается на произведение множителей вида (Ts+1), где T - вещественное число. То есть, применительно к нашим моделям, этот полином может быть представлен комбинацией интеграторов и блоков умножения. Используя при таком моделировании функцию roots, можно легко промасштабировать модель до любого порядка.
СПИСОК ЛИТЕРАТУРЫ
1. Shibata K., Yamaguchi M., Katayama H., Tokumitsu N. Mathematical Modeling for Vacuum Distillation in the Kroll Process // Nippon steel technical report. - 2002. - P. 36-39.
2. Patel M.K., Bailey C., Djambazov G., Shrimpton J., Jalili V. Application of CFD in vacuum dezincing process // Int. Conf. CFD in Minerals and Processes, CSIRO. - Melbourne, 2003. -P. 313-318.
3. Ротач В.Я. Теория автоматического управления. - М.: Изд. дом МЭИ, 2008. - 396 с.
4. Затонский А.В. Синтез систем управления сложными техническими системами // Горный информационно-аналитический бюллетень. - 2008. - № 2. - С. 82-86.
Поступила 28.02.2013 г.
UDC 66.012-52
SIMULATION OF INDUSTRIAL MANAGEMENT FACILITIES WITH DETERMINISTIC NONSTATIONARY
V.F. Bekker, N.V. Bilfeld, D.S. Kamaev*
Perm National Research Polytechnic University, Berezniki branch *Branch «Avisma» JSC «Corporation VSMPO», Berezniki
The authors have considered the task of studying the dynamics of the actual time-dependent industrial facilities by means and methods of mathematical modeling. When using standard modeling tools of Simulink library it is not possible to change the parameters of the simulated level while solving the problem. To overcome this limitation the authors proposed the method of modeling industrial objects with dynamically changing parameters in the models containing an autonomous module, which implements an arbitrary functional dependence of the parameters of the simulated standard units on machine time scale. The model is based on implementation of integrating link with a parameter that varies in time according to a certain law, or based on an external signal of any shape. The paper demonstrates the transition to the modeling of inertial element of the first, the second and higher orders and introduces the example of a two-stage relay stabilization of water temperature in a tank at various external perturbations in a flow reactor with ideal mixing. The results of computer simulation are discussed and the constraints to the model application areas are determined.
Key words:
Deterministic unsteadiness, mathematical simulations, two-position control, Simulink.
REFERENCES
1. Shibata K., Yamaguchi M., Katayama H., Tokumitsu N. Mathematical Modeling for Vacuum Distillation in the Kroll Process. Nippon steel technical report, 2002, pp. 36-39.
2. Patel M.K., Bailey C., Djambazov G., Shrimpton J., Jalili V. Application of CFD in vacuum dezincing process. Int. Conf. CFD in Minerals and Processes, CSIRO. Melbourne, 2003, pp. 313-318.
3. Rotach V.Ya. Teoriya avtomatisheskogo upravleniya [Automatic control theory]. Moscow, Publishing house MEI, 2008. 396 p.
4. Zatonsky A.V. Sintez system upravleniya slozhnymi tekhniches-kimi sistemami [Synthesis of control systems for complex technical systems]. Gorno-analitishesky bulleten - Mountain information-analytical bulletin, 2008, no. 2, pp. 82-86.
