Вычислительные технологии
Том 5, № 6, 2000
ОПТИМИЗАЦИОННЫЕ ЗАДАЧИ УПРАВЛЕНИЯ ПРОЦЕССАМИ РАЗДЕЛЕНИЯ
Н. Д. ДЕМИДЕНКО Институт вычислительного моделирования СО РАН
Красноярск, Россия e-mail: sek@icm.krasn.ru Ю. А. Терещенко Красноярский государственный технический университет, Россия
A mathematical model for the simulation of multicomponent mixtures separation processes is suggested. The problem of the optimal control has been formulated and the necessary optimality conditions have been obtained. A numerical method for the solution of optimization problems has been developed and numerical experiments have been carried out.
Введение
Исследование и проектирование химических технологий представляет собой сложную задачу, так как при ее постановке используются нелинейные системы дифференциальных уравнений в частных производных [1-6]. Математическая формулировка таких задач и вопросы их корректности, как правило, требуют специального рассмотрения. Трудности прежде всего связаны с нелинейностью уравнений и сложностью граничных условий, содержащих обыкновенные дифференциальные уравнения. С другой стороны, эти трудности обусловлены многомерностью задач, поскольку технологические процессы характеризуются довольно большим числом теплофизических и конструктивных параметров. Выбор эффективной методики решения задач моделирования и управления принято считать центральным вопросом в проблеме моделирования нестационарных режимов управляемых процессов разделения [7].
Основным инструментом проектирования технологических процессов и систем управления является вычислительный эксперимент, включающий в себя анализ многообразия технологических режимов, построение их физических и математических моделей, разработку и исследование вычислительных алгоритмов, их программную реализацию и проведение серии расчетов [8].
В настоящей работе формулируются и решаются задачи анализа математических моделей процессов разделения, задачи оптимального управления и численных экспериментов с целью проектирования энергосберегающих технологий разделения многокомпонентных смесей для промышленных объектов и создания соответствующих оптимальных систем управления.
© Н.Д. Демиденко, Ю.А. Терещенко, 2000.
1. Математическая модель управляемого процесса
Процесс разделения многокомпонентных смесей осуществляется в ректификационных колоннах на контактных устройствах (тарелках), распределенных по длине аппарата. Технологический процесс происходит в конечном числе точек объекта, однако его можно рассматривать непрерывным по длине, поэтому для моделирования возможно применение дифференциальных уравнений в частных производных. Оценка погрешности перехода от дискретной модели к непрерывной приведена в [1]. Концентрации целевого продукта в жидкости х = х(1,Ь) и паре у = у(1,Ь) — управляемые параметры, они определяются в результате решения соответствующей краевой задачи. Схема движения потоков взаимодействующих паровой и жидкой фаз приведена на рис. 1.
Разделяемая смесь в количестве Г = Г (Ь) с содержанием целевого продукта хр = хр (Ь) подается в среднюю часть колонны. В нижней ее части (кубе) происходит испарение смеси, и паровой поток V = V(1,Ь), поднимаясь вверх, контактируя со стекающей жидкостью Ь = Ь(1,Ь) и обогащаясь целевым продуктом, конденсируется в верхней части колонны (дефлегматоре) и отбирается в количестве О с концентрацией целевого продукта ха = х^(Ь). Часть сконденсированного пара Ьа = Ьа(Ь) из дефлегматора возвращается в колонну для повышения качества конечного продукта. В кубе отбирается остаток в количестве Ш = Ш (Ь) с содержанием целевого продукта хк = хк (Ь).
Одним из требований к такому промышленному объекту является его способность увеличения содержания целевого продукта в верхней части колонны и уменьшения в нижней. Важными параметрами объекта являются фазовые удерживающие способности: в колонне Нх = НХ(1,Ь), Ну = Ну(1,Ь), кубе НХк = НХк(Ь) и дефлегматоре Нхл = Нхл(¿). Индексы "х" и " у" указывают на принадлежность параметра жидкости или пару, " к" и " — кубу или дефлегматору. В колонне происходит теплообмен между жидкой и паровой средами, которые характеризуются теплосодержанием жидкости к = к(1, Ь) и пара Н = Н(I, ¿); аналогично — в кубе кк = кк(Ь), Нк = Нк(Ь) и дефлегматоре ка = ка(Ь), На = На(Ь). В куб
Рис. 1. Схема движения потоков пара и жидкости в колонне.
подводится тепло Qk, а из дефлегматора оно отводится — Qd. Коэффициент массопереда-чи ку характеризует процесс массообмена между жидкой и паровой фазами, а зависимость У* = У*(х) — равновесную концентрацию в паре. Функции х(/,£), у(/,£), х^(£) и у^)
могут быть скалярными (для бинарных смесей) или векторными (для многокомпонентных). Более подробная физическая интерпретация математической модели содержится, например, в [3].
Используя законы сохранения массы и энергии, получим математическую модель нестационарных режимов процессов разделения для колонны, куба и дефлегматора:
5(ЯхХ) д (ЬХ) = ку (У - У*) + ^ ЮФМОх,,
dt d/
d(Hyy) + _ k{v* - y) dt + d/ =ky(y y)
d (Hy H) + d(Hxh) d(Lh) + d(Vh) dV dL 5Hy dHx
----1--y +--x = Фу + Фг, (1)
d/ d/ dt dt у ^
x(/, 0) = xo(/), y(/, 0) = yo(/), V(/, 0) = Vo(/), L(/, 0) = Lo(/) (2)
с граничными условиями при / = 0
д(Hdt) = L(0,t)x(0,t) - V(0,t)y(0,t) - W(t)xfc(t), y(0, t) = (y*(xfc) - xk)a + xfc(t),
d(Hxfc hk)
dt
= L(0, t)h(0, t) - V(0, t)H(0, t) - W(t)hk(t) + Qk, dH
= L(0, t) - V(0, t) - W(t), xk(0) = Xk,o (3)
и при / = 1
dt
x)
= Vdyd - (Ld + D)xd(t), Xd(0) = Xd,o,
d (Hxd Xd)
dt
Vd(t)yd(t) - V(1,t)y(1,t) = Ld(t)xd(t) - L(1, t)x(1,t), yd(t)= y(1,t) + Ed(y*(x(1,t) - y(1,t)),
d(Hxd hd)
Vd(t)Hd - (Ld + D)hd(t) - Qd,
dHxd dt
dt
Vd(t) - (Ld(t) + D(t)), Vd - V(1,t) = Ld - L(1, t),
Vd - V (1,t)H (1,t) = Ld hd - L(1,t)h(1,t). (4)
Решение ищется в области Q = {(/, t)| / £ [0,1], t £ [0,T]}, где /, t — пространственная и временная независимые переменные соответственно; T — время управления.
Считаем в дальнейшем, что Hx = const, Hy = const, ky = k = const, Hxk = const, Hxd = const, x = x(/,t), y = y(/,t), L = L(/,t), V = V(/,t), Фн = ФН(/,t), Ф^ = Ф^(М), Фу = Фу(/,t), Ф^ = ФL(/,t) — известные функции, H = aix + a2y + a3V + a4L,
h = P\x + в2у + в?У + P4L, где ai = const, в = const, i = 1, 4. При Ed = 0 полагаем yd = y(1,t), xd = x(1,t), xd{0) = x(1, 0) = xdo, Ld = L(1,t), Vd = V(1,t).
Кроме того, Hyai + Hx^i = A, Hya2 + Hxfi2 = B, Hya3 + Hxft3 = C, Hya4 + Hxfi4 = D0. Последние допущения позволяют представить исходную систему уравнений в нормальной форме:
x't = — [Lui + xu2 + k(y - у*) + Фх(l)F(t)xF] = Xi,
Hx
xi = ui = Cl,
yt = [-VU2 - y(u4 + Фу + Фь) + k(y* - у)] = X2,
Hy
y'l = u2 =
1
VI = [АХ1 + БХ2 + п1(Уа1 - Ьрг) + щ(Уа2 - Ьр2) + п3В°+ С
+ (щ + Фу + Фь)(Н — Ьрз + Уаз) + щ(Уа4 — ВД — к) + Фн — Фн ] = Х3,
У = 4 + Фу + ФЬ = (з, Ь[ = из = Х4, Ь[ = 44 = (4. (5)
Эта система уравнений может описывать различные технологические режимы для разделения многокомпонентных смесей на промышленных объектах. Корректность некоторых таких задач рассмотрена в [1].
Система (5) решается при начальных условиях (2) и краевых условиях (3), (4). Особенность этих уравнений заключается в том, что они разрешимы относительно производных всех зависимых переменных х(1,Ь), у(1,Ь), V(¡,1), Ь(1,Ь) по независимым переменным I, Ь. Достигается это путем введения параметрических переменных , г = 1,4. Таким образом, исходная система приведена к нормальной форме. В таком виде может быть представлена любая система дифференциальных уравнений в частных производных.
Наличие параметрических переменных — основная особенность уравнений в нормальной форме. В общей записи эти переменные и управления присутствуют формально одинаковым образом, но в каждой конкретной задаче оптимального управления необходимо четко их различать, поскольку они играют принципиально разные роли: если управления могут задаваться произвольно, то значения , г = 1, 4 не задаются, а находятся по известным управлениям в результате решения задачи. Зависимые переменные непрерывны, тогда как параметрические переменные и управления в общем случае — разрывные функции независимых переменных.
Предложенная автором математическая модель используется для описания различных технологических режимов и систем оптимального управления. Корректность соответствующих краевых задач и задач оптимального управления зависит от особенностей моделируемого режима, выбора управляющих и управляемых параметров процесса и от того, полными или парциальными являются куб и дефлегматор.
2. Постановка задачи оптимального управления
В качестве управляющих воздействий из технологических соображений выбираем потоки жидкости и пара на входах в управляемый аппарат Ь(1, Ь) и V(0, Ь). На эти управляющие потоки накладываются ограничения
Ьт1п < Ь(1,1) < ЬтаХ, УтЩ < V(0, Ь) < Утах- (6)
Поскольку в дальнейшем для решения задачи оптимального управления будет использован метод вариационного исчисления, введем дополнительные управляющие функции u(t), z(t), с помощью которых ограничения (6) сводятся к равенствам
(L(1,t) - Lmin)(Lmax - L(1,t)) - U2 = 0,
(V(0,t) - Vmin)(Vmax - V(0,t)) - Z2 = 0. (7)
В качестве критерия оптимизации выбираем интеграл, характеризующий качество продуктов разделения на выходе управляемого объекта:
T 1
S = JJ (y(l,t) - 0*(l,t))2 dt ^ min (8)
0 0
(0*(/,t) — заданный состав выходного продукта).
Сформулируем следующую задачу: во множестве кусочно-непрерывных функций L(1, t), V(0, t), удовлетворяющих ограничениям (6), найти такие, которые в силу системы (1) - (4) минимизируют (8).
3. Необходимые условия оптимальности
(стационарности) в форме Лагранжа — Эйлера
Для получения необходимых условий оптимальности воспользуемся методами вариационного исчисления. Построим скалярные функции — гамильтонианы Н и Л, в области П и на границе дП соответственно:
H = (y - Г)2 + £ A.X. MiCi
i=1
i=1
(A., — множители Лагранжа);
A
(1)
h = -f- [L(0, t) - V(0, t)y(0,t) - W(t)xfc(t)] + Ak2) [y(0, t) - xfc - a(y*(xfc) - xfc)] +
H
Xfc
+Ak3)[L(0,t) - V(0, t) - W(t)] + A,
(1) V (1,t)
H
[y(1,t) - x(1,t)] + Ad2)[V(1,t) - L(1,t) - D(t)] +
Xd
+Y[(L(1,t) - Lmin)(Lmax - L(1,t)) - U2] + £[(V(0, t) - Vmin)(Vmax - V(0,t)) - Z2]
(A^, A^, Y, £ — множители Лагранжа).
Положим x = z1, y = z2, V = z3, L = z4. Рассмотрим вспомогательный функционал
J1
(y - **)2 +1 A<(X - t) +1 (z. - t
i=1 4 7 .= 1 4
d/dt
H + 2^(яГ + "яГ)
.=1
dt ö/
d/dt + / ^^ (^jdt - A.dl) z
dQ
.=1
z
Пусть Ь = а(а), I = в (а) — параметрическое задание границы д П. Тогда
31
Н + ^ (дХг +
н+ Ы -ж + -ж
г=1
<й<И + ^^ (^га'(а) — Аф'(о)) гг <а
дП
г=1
Получим вариацию 31, вызванную вариациями управлений Ь(1,Ь) и V(0,Ь):
531
п
,г=1
^ (I + д~Ж + Ж ) + £ ди. 5иг
г=1
дН
диг
¿¡¿Ь+ / ^^ (^га'(а) — Аф'(а)) 5гг <а
дп
г=1
Вариации для вспомогательного функционала 532 на границе дП вычисляются аналогично:
532
дП
дк дк дк дк дк ——-5х(0,Ь) + ——- 5у(0,Ь) + —,—т 5Ь(0,Ь) + ———- 5У (0,Ь) + — 5и+ дх(0,Ь) ; ду(0,Ь) УУ ' ; дЬ(0,Ь) ; дУ(0,Ь) 1 ' ; ди
+дк5г + ^^-5у(1,Ь) + -д^5Щ,Ь) + т дН ЛУ(1,1)
дг ду(1,Ь)
дЬ(1,Ь)
дУ (1,Ь)
<Ь+
+
дк <А
+
(1)"
5хк <Ь +
дк
+
<А
(1)
дх(1,Ь) <Ь
5х(1,Ь) <Ь — (х^Лхк + А(а1)5х(1,Ь)) |4=г.
дхк <Ь
дП х 7 дП
Таким образом, вычислена вариация вспомогательного функционала 53 = 531 + 532.
Используя аргументацию вариационного исчисления, получим следующую сопряженную задачу относительно функций Лагранжа, на основе которой разработан численный алгоритм расчета оптимальных управляющих функций [1-3]. В области П имеет место
дА1 + д^1
дЬ
д1
1I *\' I А1 I А2 , А3 ( А
к(у )Л — н + н + С н
нх Ну С \нх
Б
Ну
+
дА2 + дЬ + д1
дАз + д^з
, . А1 Аз А\ щ + Фу + Фь ща +и4 \ ~ТГ~ — +-п-а1 — в1
Нх Нх С I С С
2{у—п+к (ь — £—£ а+# Б)—А-—А- Б) —
Ну
Ну С
дЬ
д1
Нх Ну Нх С ну С
Фу + Фь (л Б \ (щ + Фу + Фь Щ., Ну Г — САз) — Аз V-С-а2 — Св2
и 2 I- Ну + Ну Б — С а2) — и4 ^С + а4) — (Фу + Фь) —
у
Ази1а1
дА4 + д^4
дЬ
д1
л и4 + Фу + Фь Щ,
с — Аз ) с аз — Свз
НТХ (А1 — Аз А) + С Р2Щ + С Щ(вз — в4) + С вз(Фу + Фь)+
, Аз _ л /и4 + Фу + Фь и* +с-и1в1 — АП -С-а4 — ~Св4
г
к
а
В этой системе неизвестные ^ или Л^ (г = 1,4) исключаются с помощью соотношений
А. Лз А V«! - ВД л + —Ь - —----Лз = о,
Н х Н х — —
^ + Н-V - - ^ - ¿02) =0,
Ну Ну — —
Л4 - ^ я
Л3
-
о,
^3 + ^4 +
А1
Нх
-ж -
Л.
Аз
-
А В
—ж - —у + Н - Л - Ь(вз + А) + V(аз + а4)
НХ Ну
Граничные условия 0 < £ < Т при I = 0 имеют вид
¿А
(1)
Л
(1)
Н
Ж(*) - Лк2) [1 + а(у*(жк)' - 1)] = 0,
Хк
Л
(1)
Л
(1)
Н
-Ь(0, *) = 0, ^2 - V(0, *) + Л
Хк
Н
(2) к
Хк
Л
(1)
gradSy = ^з - у(0, *) - Лкз) + е (К„ах + V™ - 2V(0, *)) = 0,
НХк
^4
+ тт- ж(0,£) + Лкз) = 0,
(1)
Н
Хк
2ег = 0;
при I = 1
^Л'
(1)
Л
(1)
Н
-V(1, £) - ^ (1, £) = 0,
Л
(1)
Х^
Н
-V(1, £) - ММ) = 0,
Х^
Л
(1)
Н
(у(1,£) - ж(М)) + Л^2) - Мз(1, £) = 0
Х^
gradSL = -Л^ + 7 (Ьт1п + Ьтах - 2Ь(1,£)) - ^(М) = 0,
27м = 0.
Начальные условия при £ = Т, 0 < I < 1 выглядят так:
Л1 = Л2 = Лз = Л4 = 0, Лк1) (Т) = 0, л« (Т) = 0.
Для вычисления оптимальных управляющих функций V(0,£) и Ь(1,£) применяется итерационный метод, который заключается в следующем:
1. Задаются начальные приближения V0(0,í) и Ь0(1,£).
2. Если известны ^(0, £) и Ьп(1, £), то находятся решения прямой и сопряженной задач.
3. Полагаем Vra+1 (0,*) = Vra(0,í) - TlgradSy, Ьга+1(1,£) = Ьп(1,£) - T2gradSL.
4. Предельные значения Ь(1,£) и V(0,£) дают решение задачи оптимального управления.
У
к
к
к
0
4
Х^, Хк, мол. доли 2 „_
0.8 -- 1/
0.6 - /
0.4 -1-1-'-'-
0 5 10 15 20 ч
Рис. 2. Концентрация бутана в дефлегматоре в пусковом режиме при управлении Ь(1,Ь). Кривые 1, 2 — начальное и оптимальное управление соответственно.
Для численного решения краевых задач и задач оптимального управления разработан численный алгоритм с использованием метода центральных разностей и треугольных сеток [1].
В качестве примера приведены результаты расчетов оптимальной управляющей функции при оптимизации пускового режима для промышленной колонны К-34 установки сернокислотного алкилирования изобутана бутиленами (разделяемая многокомпонентная смесь сведена к бинарной). На рис. 2 показано изменение концентраций целевого продукта в дефлегматоре в переходном режиме, а на рис. 3 изображена оптимальная управляющая функция ¿(1, ¿). При оптимальном управлении выход на заданное значение концентрации целевого продукта происходит быстрее, чем при неоптимальном. Подробное описание установки и экспериментальные значения основных параметров процесса приведены в [3].
Заключение
Разработан метод математического моделирования нестационарных режимов разделения многокомпонентных смесей для исследования и проектирования систем оптимального управления ректификационными колоннами. Метод апробирован на промышленных ректификационных установках. Развитая в работе общая теория и метод анализа нестационарных режимов могут быть применены к широкому классу технологических аппаратов: колоннам ректификации (насадочным и тарельчатым), абсорберам, теплообменникам и др. Предлагаемый метод анализа динамических характеристик объектов с распределенными параметрами прошел экспериментальную проверку. С помощью предложенных вычислительных алгоритмов исследованы возможности оптимизации пусковых режимов, выполнено моделирование перехода от одного стационарного режима работы к другому, стабилизации заданного состава выходных продуктов.
Список литературы
[1] ДЕМИДЕНКО Н. Д. Моделирование и оптимизация тепломассообменных процессов в химической технологии. М.: Наука, 1991.
[2] ДЕМИДЕНКО Н. Д. Управляемые распределенные системы. Новосибирск: Наука, 1999.
¿(1,4), кмоль/ч
50 0
10 15 20 ч
Рис. 3. Оптимальная управляющая функция Ь(1,Ь) в пусковом режиме.
5
[3] Демиденко Н. Д., Ушлтинскля Н. П. Моделирование, распределенный контроль и управление процессами ректификации. Новосибирск: Наука, 1978.
[4] DEMIDENKo N. D. Modelling of optimal regimes in chemical engineering objects with interacting flow recirculation // Syst. Anal. Model. Simul. 1987. V. 4. P. 309-320.
[5] DEMIDENKo N. D. Optimal control of the complicated objects distributed parametres // Syst. Anal. and Simul. 1985. V. 27, No. 1. P. 425-428.
[6] DEMIDENKo N. D. Problems on optimisation of information measuring systems with distributed parametrs // Syst. Anal. Model. Simul. 1990. V. 11-12. P. 907-920.
[7] Кафаров В. В., ВЕтохин В. Н. Основы построения операционных систем в химической технологии. М.: Наука, 1980.
[8] Вычислительный эксперимент в проблеме цунами / Ю. И. Шокин, Л. Б. Чубаров, Ан. Г. Марчук, К. В Симонов. Новосибирск: Наука, 1989.
Поступила в редакцию 4 мая 2000 г., в переработанном виде — 20 июня 2000 г.