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

УДК 681.5.01

Н. Д. Демиденко, Ю. А. Терещенко, И. Н. Мельник

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ И ОПТИМИЗАЦИЯ СИСТЕМ С РАСПРЕДЕЛЕННЫМИ ПАРАМЕТРАМИ

Решается задача оптимального управления процессами разделения многокомпонентных смесей. Получены необходимые условия оптимальности. Разработан численный метод решения задачи оптимизации. Проведен численный эксперимент.

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

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

Уравнение нестационарного процесса массообмена для многокомпонентных смесей. При исследовании процесса ректификации интерес представляет распределение концентраций компонентов по длине колонны в статических и динамических режимах ее работы. Поэтому математическая модель процесса ректификации, как правило, представляет собой систему уравнений, записанных относительно концентраций компонентов. Примером такой математической модели является следующая система уравнений [1]:

д(Нхх') д(Ьх') = к (у _/)+ 0

д д1 № У') Р!^

+даТ=ку(_ У')+р Ч’ 1 - ' - *’ (1)

где ж, = х' ( (), у' = у' ( t) - концентрации компонентов в жидкой и паровой фазах; Нх, Ну - удерживающие способности по жидкости и пару; Ь, V - потоки жидкости и

* * / \ ЧУ

пара; у = у (х) - равновесная концентрация в жидкой фазе; ку - коэффициент массопередачи.

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

Влияние гидродинамических явлений на процесс мас-сообмена можно выразить зависимостью коэффициента массообмена от величины парового потока ку = кУ, где к определяется по экспериментальным данным [1].

Такая зависимость является приемлемым приближением для некоторых режимов работы ректификационной колонны. Равновесная концентрация в паре определяется по концентрацию компонентов X; эмпирической зависимостью у* = ( Т2 + Р2Тк + Рз'), где Р - давление в колонне; Тк - температура кипения жидкости; Ри, Р2,, Рм - коэффициенты, полученные по эмпирической зависимости давления чистых компонентов от температуры методом наименьших квадратов, ' = 1, N.

Рассмотрим ввод потоков сырья в ректификационную колонну не по всей длине, а в некоторой ее части, которую будем называть областью ввода. Функции р1 / (I, t),р2^ (I, t), определяющие плотности потоков г-го компонента по длине колонны, находятся следующим образом:

р1Л = Рь (t()х)\(1/1 ()’р2л = ру(t)уf (1/2(), гдеРЬ (t)’Ру ^) - потоки сырья в жидкой и паровой фазах; х^ ^), у^ ^) - концентрация г-го компонента в сырье. Система (1), как было показано в [1], может быть приведена к виду

дРи _ с ЭРи дt 1 у'

Эр2

■+ С-

- = ку ( - У*)+ Рі — = ку - У,)+ Р 2/1, 1 < і < N■

(2)

дґ '2 ' У'

Таким образом, имеем систему уравнений в частных производных с постоянными коэффициентами. Функции хі и Уі определяются теперь следующим образом: х = Ріі у = Р 2і і N 5 Уі N ■

X Ри X Р2і

і=1 і=1

Система уравнений (2) описывает процесс массооб-мена, происходящий внутри ректификационной колонны. Кроме этого, необходимо учесть процессы, протекающие в кубе и дефлегматоре, которые находятся в концах колонны. Уравнения, описывающие процессы в кубе и дефлегматоре, определяют граничные условия для системы уравнений (2).

Начальные условия имеют вид

Ріі (І, 0) = Фіі (І), р 2і (І, 0) = Ф2і (І), 1 < і < N . (3)

Постановка краевой задачи. В большинстве случаев процесс ректификации осуществляется с рециркуляцией выходных потоков в кубе и дефлегматоре. Граничные условия для этих уравнений задаются на границе І = Ь , І = 0. Значения концентрации компонентов в дефлегматоре определяются по уравнению покомпонентного материального баланса для дефлегматора

* (Н*хЛ (0),

Л

■ = С2Р2і (Ь 0 - (0х*і ()- П(І)х*і ()■ (4)

Начальные условия для уравнения (4) следующие:

Фц (ь)

■(0) = -

£фц(ь)

(5)

1=1

Аналогично задаются граничные уравнения для остальных N уравнений системы (2). Граничные условия задаются также в точке / = 0.

Концентрации компонентов жидкости, находящейся в кубе, определяются по уравнениям покомпонентного материального баланса для куба:

Ж Я*,х*(/}) = с.р,(0, /) - ж(/(/)_ у(/)(о, /) (6)

Фк(о)

(7)

п+1т пт п+1т+1 п+1т

р1/ - р1/ __ с рЦ - рЦ = кпт1

Д/

Д/

п+1т ¡-.пт п+1т+1 п+1т

р2/ р2/ , с р2/ р2/ = -пт (_^*пт упт) + р пт

~ 1 К1 ~ У V/ У/ р2

Д/

д/

к7 = кс2£р2т, уТт = (р1/ )2 + Р2Т"т + Р3/ \"т1рт ,

¡=1

Т^т в пт + ^(в пт )2 ^2дпт

д„т = £р х™ впт = £р .х™ С™ = £Р2]х"т _ Рт ¡=1 1 1 ¡=1 ¡_1

Рт = Ря _(я _ Рв )тМ/Ь,

/Ы / N

пт пт пт пт

£р1 , У/ = Р2у £Р2| ,

где £ - длина объекта.

Начальные условия при / = 0 примут вид р™ = ф|”, рпт = фт , граничные условия при / = 0 перепишутся следующим образом:

р2+11 = у^у^К, уГ = (( + Ратя1 + Рз,)+1 /Р1,

1 хк хк я хкхк _ Л _«1 т/и,,и тттп_.п , -V.1 = Ф1/

17«+^«+1 _и п V«

п1 п п п п ,

= с1р1/ _ У0ук _ Ж хы’

£фи-

1=1

УЦ = С1 £ рп _

при I = Ь:

¡=1

и+Ш ги+]

Р1/ = Ь

7 с1 >

£7 м+^И+1 _ и п уп N

11 хЖ лсИ 11 хЖл<й _ _ ппМ „п ^ иМ = с2р 2/ ХЖ/с2 £ Р

фМ х1 = _%_

ХЖ/ N

ЯМ+1 Ц!

хЖ ЯхЖ

2 ^ р 2| 5 ¡=1

£ф]

¡=1

А/

\ ' _мМ тп пп 7^1 = С2 £ р 2| _ ЬЖ _ О , ЯхЖ = с3-¡=1

£ф^^ (0)

¡=1

Итак, сформулированная краевая задача (2).. .(7) является математической моделью работы ректификационной колонны без системы регулирования в нестационарном режиме. Решение этой краевой задачи не всегда имеет технологический смысл, который для ректификационной колонны соответствует аварийному режиму работы. При этом происходит залив колонны или ее исчерпывание. С помощью системы регулирования поддерживаются такие значения входных и выходных потоков, чтобы выполнялся внешний материальный баланс для колонны. В математической модели это можно отразить следующим образом: ^ (/) + ЕУ (/) = О(/) + ж (/).

Численный метод решения краевой задачи. Так как решить аналитически сформулированную краевую задачу не удается, то для этого следует применить численный метод [2]. Рассмотрим алгоритм численного решения краевой задачи (2).. .(10). Для этого область решения покроем равномерной сеткой с шагом Д/ по временной координате и с шагом Д1 по пространственной. Для решения уравнений на сетке выбран трехточечный шаблон.

Заменим производные в уравнениях соответствующими конечно-разностными соотношениями:

пт пт *пт пт

- = ку ( _ у,- )+ р^ >

Для определения коэффициента к используются экспериментальные данные. В работе [3] приведены значения концентраций в выходных продуктах для двенадцати статических режимов колонны К-34. Для каждого режима найдено такое значение коэффициента к, при котором расчетные значения выходных концентраций имеют наименьшее отклонение от экспериментальных данных, и взято среднее арифметическое значение этого коэффициента. Среднее значение коэффициента для К-34 равно 0,395, что близко к проектным значениям КПД тарелок. Разброс значений для коэффициента можно объяснить следующими факторами:

- недостаточной точностью измерения концентрации;

- получением экспериментальных данных не в статическом режиме (для некоторых режимов не соблюдается покомпонентный материальный баланс);

- неполной адекватностью математической модели реальному процессу.

Оптимальное управление в замкнутой системе с измерением параметров в точках, распределенных по длине аппарата. Применить на практике непрерывное измерение удается в лишь редких случаях, поэтому приходится пользоваться измерением параметров в некоторых точках. В связи с этим возникает вопрос нахождения оптимальных координат точек измерения и весовых коэффициентов. Рассмотрим решение такой задачи при заданном числе точек измерения. Запишем для случая дискретного измерения уравнение обратной связи [4; 5]

I(/) = «°(/)'+ ££ ■КI а к ( г >

к=1 5=1 0

_с( чк)2

с11.

(8)

где и! (/),и! (/) - управляющие функции; R - количество точек измерения к-й измеряемой функции j-й управляющей функции; - координата ^-й точки к-й измеряемой функции j-й управляющей функции; к к - весовые коэффициенты; а к (/,/) - контролируемый параметр.

На координаты точек измерения накладываются ограничения:

0 </к < Ь, 1 < j < 8, 1 < к < 3, 1 < 5 < К . (9)

Задачу оптимального управления в замкнутой системе с дискретным измерением можно сформулировать следующим образом: во множестве действительных чисел необходимо найти такие значения параметров кк и /к, при которых решение системы уравнений (2)...(7) дает минимум функционалу 5 с учетом ограничений (9):

5 = £1[ (/Ь ви )2 + к2/ (хй(/) _ 6 2/ )2 ], (10)

¡=1 0

где к1/, къ - коэффициенты, определяющие ценность г-го компонента; 91(, 9 2/ - заданные числовые значения концентраций выходных продуктов; хы (/), хл (/) - состав выходных потоков; Т - время управления.

х

Ж

и

Необходимые условия оптимальности содержат следующую сопряженную задачу:

_ с дк=її. дґ 1 ді РЬ

( _п,)р,+£А +£(?і _л¡)Рі

р,Ь _ Рсі

' Ь

(11)

дґ

ді

_ І )+Ёіі _ Лі)

с начальными условиями

/ 8 Л

І і (і, т) = «2і (і, Т )££ю і (Т) (і, ґ) = 0

1=1 ¿=1

и граничными условиями при і = 0,0 < ґ < Т:

(12)

0,ґ) = £ Хкі(ґ)

і=1

і 4,(0

нгі

^ а3і (0, ґ)£)г(0)

сг

1=1 г=1 сі

>4і(ґ) = с2Лі (0,ґ),

АЙ (ґ) , 1 ^ Ь (ґ) _ М1 (ґ) _ М4 (ґ) + Му (ґ) +

------------" = ~к к (ґ)------------------------------------------------------------+

Іґ

+ & ( ) 0 ( ) + 2к2і [[ (ґ) _ ^2і (ґ)]] (Т) = 0; Нхк

2

(13)

при і = Ь, 0 < ґ < Т :

Лі ґ) = -

і=1

Ні

8 Л

- + ^ (ґ)діі (ґ) =

= с1 §і (Ь, ґ) _ ££ ші (ґ)«зі (Ь, (Ь),

1=1 Г=1

^ (ґ)м8(ґ) + + 2кц [хц (ґ) _ 01і (ґ)]іі (Т) = 0,

аґ

Н

ме с наиболее эффективной управляющей функцией Б(ґ) (величина отбора верхнего продукта) для дискретного измерения параметров.

На основе разработанных алгоритмов проведено численное исследование систем оптимального регулирования при различных типах регуляторов и различных возмущениях. В качестве объекта управления взята колонна К-34 (дебутанизатор) установки сернокислотного алки-лирования изобутана бутиленами. Расчеты проводились для двух видов возмущений по составу сырья у!2 (рис.1); ступенчатого (рис. 1, а) и линейного (рис. 1, б).

Х/ 2 0,6 %-

0,4

0,2

_1_

_1_

_1_

5 10 15 20 ґ, Ч. а

0,5

0,3

_1_

_1_

_1_

5 10 15 20 ґ, Ч.

^ ^) с2р 2і (Ь,ґ)2 Хіі(ґ)Уо(ґ),13і (т ) = 0,

іґ

і=1

НХ

(14)

где

А „ = и і (ґ)

і (і) _£ і (і)

к=1 ІІ дР1і (і,ґ) іі ді 1

ні

д& і (і, ґ)а 2і (і, ґ)

дґ

к _с(і_к)2

4 (/) = к>

Таким образом, задача оптимизации управления в замкнутой системе с дискретным измерением сводится к решению системы уравнений (2).. .(8), (11)_(14) с уче-

том (9). Эта задача решается численным методом итерационно. При некоторых начальных значениях парамет-ров/к]г икк]г две краевые задачи (2)...(7)и (11)...(14) решаются последовательно методом конечных разностей. Затем находятся значения величин ) и ) . Новые зна-

^ ¿г,к ' 1 '/

чения параметров | и к1р определяются по формулам

■г )/+1 = кг У _ я Г« )к', (/- У+1 = (/- У _ я кГ 2)(<)/',

где г - номер итерации; , Л|.2) - коэффициенты, которые выбираются такими, чтобы на каждой итерации происходило уменьшение функционала (10).

Численное решение задачи оптимизации. Приведем в качестве примера решение задачи стабилизации состава выходных продуктов при возмущении по составу сырья. Суть задачи следующая: при работе ректификационной колонны в статическом режиме возникает возмущение по составу сырья и необходимо с помощью системы управления поддерживать постоянный состав выходных потоков, такой же, как в статическом режиме. Рассмотрим решение задачи при управлении в замкнутой систе-

Рис. 1. Возмущения по концентрации бутана в сырье

Начальные значения коэффициентов усиления регуляторов следующие: для П-регулятора вверху колонны К0 = 9 000, внизу К0 = 15 000; для И-регулятора вверху колонны К0 = 0,05, внизу К0 = 0,05 ; для Д-регулятора вверху колонны К0 = 28 000 , внизу К 0 = 28 000 . Оптимизация процесса управления осуществляется с помощью двух регуляторов внизу и вверху колонны. Соответ-

-I-----------1---------1---------1--------1---------1--------1---------1--------1---------1--------1---------1—

"В*

2 4 Б

10 12 I, час

Рис. 2. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с П-регулятором: а - вверху колонны (х° и х^О); б - внизу колонны (х° и х^пт(/)}; в - управляющего потока (О0 и Оопт(/))

і=1

і=1

х

ха

І=1

ственно оптимизируемый функционал содержит квадраты отклонений концентрации компонентов вверху и внизу колонны.

Приведем начальное и оптимальное изменения управляющего потока отбора верхнего продукта О0,Оопт (/) в системе управления с П-регулятором (рис. 2). Начальное значение показателя качества одинаково для всех регуляторов при ступенчатом возмущении и равно р = 17,334 76 = 0,77219 +16,562 57 . Первое слагаемое в последней сумме означает часть показателя качества для верхнего продукта, второе слагаемое - для нижнего продукта. Оптимальные значения коэффициента усиления для П-регулятора К = 10 500 , К ™т = 18 000, ропг = 0,004 394 = 0,004 310 + 0,000 084.

Оптимальные и начальные значения управляемых и управляющего параметров в системе с И-регулято-ром (рис. 3) следующие: К™т = 1,90,К™т = 1,90 ,

ропт = 0,679301 = 0,014596 + 0,664705. Отметим, что у системы с И-регулятором показатели качества хуже, чем у системы с П-регулятором.

а *0

■ПО

4 6 б 8 10 12

—I—I—I—I—I—I—I—I—I—I—I—I-------->

О 2 4 6 В 10 12 1, час

Рис. 3. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с И-регулятором: а - вверху колонны (х° и х“1 ; б - внизу колонны (х° и х“ 11 (/));

в - управляющего потока (О0 и О опт(/))

Представим результаты расчетов системы регулирования с Д-регулятором (рис. 4):

ропт = 0,013482 = 0,012405 + 0,001077, К0™ = К= 329000.

Получены следующие результаты расчетов оптимальных параметров с ПИ-регулятором (рис. 5): К 0Пт = 10500, КОТ = 18000, К0Пт = 0,93, К0Пт = 0,91 , ропт = 0,004375 =

= 0,004244 + 0,000131. Управление в этом случае более качественное, чем в системах отдельными П- и И-регулято-рами.

Рис. 4. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с Д-регулятором: а - вверху колонны (х°л и х™ (/)); б - внизу колонны (х°к и х™т(/)); в - управляющего потока (О0 и Оопт (/))

Рис. 5. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с ПИ-регулятором: а - вверху колонны ( х0 и х°Л (/)); б - внизу колонны (х0 и хГ (/)); в - управляющего потока (О0 и Оопт (/))

Приведем результаты расчетов с ИД-регулятором (рис. 6):

ропт = 0,02294 = 0,00888 + 0,014061, К0™ = К0™ = 331000.

..

0,87'

0,86.

0,85'

0,84.

0,83,

Ч-----------1--------1--------1---------1-------1---------1--------1-

Рис. 6. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с ИД-регулятором: а - вверху колонны (х° и х°/Т (/)); б - внизу колонны (х0 и х™т(/)); в - управляющего потока (О0 и Оопт (/))

а

50 ■

40 ' 30 .

20 ■

10 ■

!У" I'/

-I—I—I—I-

-I—I—I—I-

•—1------->

О 2 4 6 в ю 12 1, час

Рис. 7. Начальные и оптимальные значения концентраций компонентов целевого продукта при ступенчатом возмущении в системе с ИД-регулятором: а - вверху колонны (х° и х^ш(/)); б - внизу колонны (х° и хЦш(/)); в - управляющего потока (О0 и О опт(/))

Результаты расчетов с ПИД-регулятором следующие (рис. 7): ропт = 0,004378 = 0,004297 + 0,000081. К ОЦт = 10500,

К™т = 18000, К™т = 0,03,К™т = 0,01, К°™ = 75000, К™т = 25000 . Регуляторы П, ПИ, ПИД равнозначны по результату регулирования при управлении качеством верхнего и нижнего продуктов.

Приведем результаты расчетов систем авторегулирования при линейном возмущении по составу сырья (см. рис. 1, б) с регуляторами П, Д, ПИ: р0 = 5,485232 =

; = 0,178581 + 5,306651. Для П-регулятора ропт = 0,009916 = = 0,001835 + 0,008081, Ко™ = 84000, К^Т = 28000. Для Д-регу-лятора ропт = 0,119634 = 0,006448 + 0,113186 , К0™ = 240000,. К<отп = 250000 (рис. 9). Для ПИ-регулятора ропт = 0,001856 =

о К п

= 0,2, К по

0,05 .

к

0,5

0,4

0,3

0,2

0,1

Ч------------1---------1---------1---------1---------1---------1---------1----------1--------1----------1-

10

12 час

Рис. 8. Начальные и оптимальные значения концентраций компонентов целевого продукта для систем авторегулирования с П-регулятором: а - вверху колонны ( х0 и х™ (/)); б - внизу колонны ( х0 и хГ (/)); в - управляющего потока (О0 и О опт(/))

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

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

= 0,001383 + 0,000472, К о„ = 14000, К опт = 18000,

ния оптимального управления. Число итераций не превышает трех.

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

1. Демиденко, Н. Д. Управляемые распределенные процессы / Н. Д. Демиденко. Новосибирск: Наука. Сиб. отд-ние, 1999. 394 с.

2. Демиденко, Н. Д. Оптимизационные задачи управления процессами разделения / Н. Д. Демиденко, Ю. А. Терещенко // Вычислительные технологии. 2000. Т. 5. № 6. С. 36-44.

3. Демиденко, Н. Д. Моделирование, распределенный контроль и управление процессами ректификации / Н. Д. Демиденко, Н. П. Ушатинская. Новосибирск: Наука. Сиб. отд-ние, 1978. 288 с.

4. Demidenko, N. D. Modelling jf optimal regimes in chemical engineering objects with interacting flow recirculation / N. D. Demidenko // Syst. Anal. Model. Simul. 1987. Vol. 4. P. 309-320.

5. Demidenko, N. D. Problems on optimization of informationmeasuring systems with distributed parameters / N. D. Demidenko // Syst. Anal. Model. Simul. 1990. Vol. 11-12. P. 907-920.

t. час

Рис. 9. Начальные и оптимальные значения концентраций компонентов целевого продукта для систем авторегулирования с Д-регулятором: а - вверху колонны (х° и х^1™ ()); б - внизу колонны (х° и х“т ()); в - управляющего потока (О0 и О опт ())

Рис. 10. Начальные и оптимальные значения концентраций компонентов целевого продукта для систем авторегулирования с ПИ-регулятором: а - вверху колонны (х° и х^1”” ()); б - внизу колонны (х°° и х^шг(/)); в - управляющего потока (О0 и Оопт ())

N. D. Demidenko, J. A. Tereschenko, E. N. Melnik

MATHEMATICAL MODELING AND OPTIMIZATION OF SISTEMS WITH DISTRIBUTED PARAMETERS

It is stated the problem of optimal control processes of multicomponent mixtures division. It is obtained the necessary optimal conditions. It is developed the numerical method of optimization problems solving. It is carried out numerical experiment.