Численный алгоритм поиска оптимального температурного режима химического процесса на основе эволюционных вычислений Текст научной статьи по специальности «Математика»

ДИФФЕРЕНЦИАЛЬНЫЕ УРАВНЕНИЯ И

ПРОЦЕССЫ УПРАВЛЕНИЯ N. 2, 2024 Электронный журнал, рег. Эл № ФС77-39410 от 15.04.2010 ISSN 1817-2172

http://diffjournal. spbu. ru/ e-mail: _ jodiff@mail.ru

Компьютерное моделирование динамических и управляемых систем

Численный алгоритм поиска оптимального температурного режима химического процесса на основе эволюционных вычислений

Антипина Е.В.*, Мустафина С.А.**, Антипин А.Ф.***

Уфимский университет науки и технологий

*

stepashinaev@ya.ru

**

mustafina_sa@mail.ru

***

andrejantipin@ya.ru

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

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

1. Введение

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

позволяет исследовать закономерности протекания процессов без проведения лабораторных экспериментов.

Одной из таких задач является задача оптимального управления химическим процессом, постановка которой включает в себя три основные части [1].

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

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

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

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

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

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

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

2. Постановка задачи

Пусть химический процесс описывается системой обыкновенных дифференциальных уравнений

йх

— = Г(х(1),Т,1) (1)

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

х(0) = х0, (2)

где х(Ь) = (х1(1),х2(1),...,хп(1))т - вектор концентраций веществ, Т = Т(Ь) - температура реакционной смеси, t Е [0, т] - время,

/(х(£),Т,1) = (/1(х(£),Т,1),/2(х(£),Т,1),..., /п(х(1),Т,1))т - вектор-функция, непрерывная вместе со своими частными производными.

В системе дифференциальных уравнений (1) с начальными условиями (2) переменными состояния являются концентрации веществ х^) (I = 1,п), управлением - температура реакции Т.

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

т<т(г)<т,гЕ[0,т], (3)

где Т_, Т - нижняя и верхняя допустимые границы температуры, которые устанавливаются ввиду технологических ограничений.

Требуется определить оптимальный температурный режим Т(1) химического процесса, описываемого системой дифференциальных уравнений (1) с начальными условиями (2) с учетом ограничений (3) доставляющий максимум функционалу

/СО = (р(х(т)). (4)

Функционал (4) может выражать выход целевых продуктов в конце реакции, содержание побочных веществ, степень превращения реагентов, селективность и другие показатели эффективности химического процесса.

Для построения аппроксимирующей задачи на интервале [0, т] введем сетку дискретизации с шагом Ь = — и узлами ...^н, такими, что 10 < < ■■■ < , 10 = 0, = т, в которых

будем искать управляющую функцию Г(0. Для получения промежуточных значений управления используем кусочно-постоянную аппроксимацию

Т(г)=Т(1к) = Тк, 1Е[1к,1к+1], (5)

где к = 0, N - 1.

Ограничение (3) преобразуется к виду

Т<Тк<Т, к = 0,М-1. (6)

Конечномерная задача, аппроксимирующая задачу (1)-(4) заключается в определении температурной последовательности (5), при котором с учетом ограничений (6) функционал

КТ) = у(х(1м)) (7)

Дифференциальные уравнения и процессы управления, N. 2, 2024 достигает наибольшего значения.

3. Алгоритм расчета оптимального температурного режима химического процесса

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

Пусть математическим аналогом особи является вектор температуры Т = (Т0,Т1,...,Ты-1), где Тк = Т^к), к = 0, N — 1. Популяцией является набор из Р векторов Т = (Т]0, Т^1,.. .,Т]н-1),

] = тр.

Для поиска наилучшего решения используется значение функции приспособленности, которое показывает, насколько хорошо подходит особь в качестве решения задачи. В задаче поиска оптимального температурного режима в качестве функции приспособленности будет выступать функционал (7). Чтобы вычислить его значение, необходимо найти численное решение системы дифференциальных уравнений (1) с начальными условиями (2).

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

Шаг 1. Сгенерировать начальную популяцию векторов температуры Т° = (Т^, Тд,..., Т°ы-1),

] = ЪР:

Т°к:=Т + Ч(Т — Т),

где ц Е [0,1] - случайное число.

В качестве текущей популяции задать начальную популяцию: 1 = 0.

Шаг 2. Вычислить значение функционала (7) для каждого вектора температуры Т--,} = 1, Р. Шаг 3. В качестве вектора-мишени задать первый вектор текущей популяции: Тт = Т-. Шаг 4. В текущей популяции найти вектор температуры Ттах, которому соответствует наибольшее значение функционала (7).

Шаг 5. Выбрать случайным образом четыре различных вектора Та, Ть, Тс, Та, которые не являются вектором-мишенью Тт и вектором Ттах Создать новый вектор-мутант по правилу:

ТтШ к' = Ттах к + ^(Так + ТЬк — Тск — Тйк),

где г Е [0,5; 1] - параметр мутации [13], к = 0^ — 1. Если Ттш к < Т или Ттш к > Т, то Ттш к присвоить случайное значение температуры из отрезка [Т, Т].

Шаг 6. Создать вектор-образец ТоЬг в результате применения операции кроссовера к вектору-мишени и вектору-мутанту. Для этого сгенерировать для каждой к-й координаты (к = 0, N — 1) случайное число рк из отрезка [0; 1]. Координаты вектора-образца определяются по формуле:

Т _ (ТтШ к, ^ 9, 1оЬГк = \Ттк, Кк>9,

где д Е [0; 1] - параметр скрещивания [13].

Шаг 7. Вычислить значение целевого функционала (7) для вектора-образца ТоЬг. Если /(ТоЬг) > /(1т)), то в новую популяцию вместо вектора-мишени поместить пробный вектор, т.е. Тт = ТоЬг. Иначе вектор-мишень Тт остается в популяции.

Шаг 8. Если в качестве вектора-мишени рассмотрены все векторы температуры текущей популяции (Тт = Тр), то перейти к шагу 9. Иначе Тт = Т-1п+1 и перейти к шагу 4.

Шаг 9. Проверить условие окончания поиска решения. Вычислить расстояние р^ (I,} = 1, Р) между векторами температуры текущей 1-й и предыдущей (I — 1)-й популяций, а также изменение значений целевого функционала А^ (I > 1):

N-1

PiJ = р(Т1,Т}-1) = \\Tl - tJ-1\\ = - Tj-1)2

к=0

Atj = AU(Tll),J(TJ1-1)) = \КТг1) -J(Tf-1)\, i,j = 1,P.

Если на протяжении г поколений выполнены условия р^ < г, А^ < г, то есть происходит незначительное изменение популяции и целевого функционала, то необходимо остановить поиск решения и выбрать из последней популяции вектор температуры' которому соответствует наибольшее значение целевого функционала (7). В противном случае увеличить номер популяции I на 1 и перейти к шагу 3.

4. Вычислительный эксперимент

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

Х- ^ Х2, Х2 ^ Х4, Х- ^ Х3, Х- ^ Х4, Х2 ^ Х3, Х3 ^ Х$, (8)

где Х1 - нафталин, Х2 - нафтохинон, Х3 - фталевый ангидрид, Х4 - углекислый газ, Х5 -малеиновый ангидрид.

Кинетические уравнения согласно закону действующих масс имеют вид:

Т) = к1(Т)х1, ^2 (%, Т*) = к2(Т)Х2, ж3(х,Т) = к3(Т)х1, ж4(х,Т) = к4(Т)х1, ж5(х,Т) = к$(Т)Х2, ж6(х,Т) = к6(Т)Хз,

где №](х,Т} - скорость у'-й стадии реакции (8) (моль/(л-ч)), } = 1,6; хг - концентрация /-го вещества (моль/л), I = 1,5; - константа скоростиу'-й стадии (1/ч), } = 1,6, зависящая от

температуры Т исходя из уравнения Аррениуса

kj(T) = kojexp(-^),

где к0у - предэкспоненциальный множитель (1/ч), Еу - значение энергии активации у'-й стадии (Дж/моль), Я - универсальная газовая постоянная (8.31 Дж/(мольК)).

Численные значения кинетических параметров реакции получения фталевого ангидрида приведены в работе [14].

Динамику концентраций веществ реакции получения фталевого ангидрида можно представить в виде системы обыкновенных дифференциальных уравнений:

= —]к1(х, Т) — Wз(x, Т) — ]м4(х, Т), = W1(x, Т) — W2(x, Т) — ws(x, Т),

йх3 М

-¿3 = w3(x, T) + ws(x, T) - w6(x, T), (9)

^4 = w2(x, t) + w4(x, T),

с начальными условиями xt(0) = xf.i = 15.

Пусть на значения температуры наложены технологические ограничения: 620К<Т< 644К (10)

Продолжительность процесса примем за т = 0.8 ч. Начальные концентрации веществ заданы значениями (моль/л):

х^0) = 1,х^(0) = 0.i =25. (11)

В качестве критериев оптимальности можно выбрать показатели эффективности проведения процесса на стадии его завершения (t = т):

1) максимальная степень превращения исходного вещества - нафталина:

(1 — х1(т)) • 100% ^ max. (12)

2) максимальная концентрация целевого продукта реакции - фталевого ангидрида:

х3(т) ^ max. (13)

3) минимальная концентрация побочного продукта - углекислого газа:

—х4(т) ^ max. (14)

Для решения поставленных задач создана программа в среде визуального программирования Delphi. Разработанный алгоритм дифференциальной эволюции применен со следующими параметрами: размер популяции Р = 70, параметр мутации z = 0,7, параметр скрещивания д = 0,8, параметры завершения поиска решения г = 5, г = 10-5. Количество точек разбиения временного интервала N принято равным 100. Численное решение системы дифференциальных уравнений (9) с начальными условиями (11) получено с помощью явного метода Рунге-Кутта четвертого порядка точности.

В результате расчетов установлено, что для достижения наибольшей конверсии нафталина, равной 99,9%, необходимо вести процесс при изотермических условиях, поддерживая температуру реакции, равной верхней допустимой границе (644 К). На рис. 1 приведено изменение степени превращения исходного вещества во времени при найденном температурном режиме.

0 4---1-■-1-.-1---1—

0,0 0,2 0,4 0,6 0,8

г, ч

Рис. 1. Динамика степени превращения нафтохинона

Для получения наибольшей концентрации фталевого ангидрида следует начать процесс при температуре, равной 644 К, и поддерживать ее на этом уровне на протяжении 0,1 ч. Затем до конца реакции необходимо постепенно снижать температуру до 622 К (рис. 2). При этом максимальная концентрация целевого продукта реакции составляет 0,73 моль/л (рис. 3).

620-1-■-1-■-I-■-т-■-1—

0,0 0,2 0,4 0,6 0,8

/, ч

Рис. 2. Субоптимальный температурный режим

0,8 и

0,6-

"л

I °'4"

А

К

0,20,00,

Рис. 3. Динамика концентрации фталевого ангидрида

Чтобы обеспечить наименьшее содержание побочного продукта реакции, нужно начать процесс при минимальной допустимой температуре (620 К) и вести процесс при данной температуре до его завершения. При данном температурном режиме наименьшее значение концентрации углекислого газа составляет 0,18 моль/л (рис. 4).

0,20-,

0,15-

EJ Л

г; о

0,10-

0,05-

0,00-1 0,

Рис. 4. Динамика концентрации углекислого газа

Численные решения задач поиска оптимального температурного режима с целевыми функционалами (12), (13), (14) найдены также с помощью метода вариаций в пространстве управления [15], реализованного авторами в программе на языке программирования Delphi. Задача решена с шагом по управлению, равным h = 0,01, и количеством точек разбиения интервала времени N = 100.

Решение задачи с целевым функционалом (12) рассчитано при двух начальных приближениях управления, равных 632 К (середина интервала допустимой температуры) и 620 К. В обоих случаях получен аналогичный изотермический режим (644 К), как режим, вычисленный с помощью алгоритма дифференциальной эволюции.

Для решения задачи оптимального управления с целевым функционалом (13) использовалось три начальных приближения: 622 К, 634 К и 644 К. Результаты вычислений показали, что при

Т-1-1-'-1-'-г

0 0,2 0,4 0,6 0,8

0 0,2 0,4 0,6 0,8

U ч

начальных значениях температурной зависимости, равных 622 К и 634 К, температурная кривая по структуре совпадает с температурным профилем, вычисленным с помощью алгоритма дифференциальной эволюции. При начальном приближении управления, равном 644 К, решение попало в локальный экстремум. Наибольшая концентрация фталевого ангидрида составила 0,71 моль/л при соблюдении температурного режима, приведенного на рис. 5.

644 640 636 ^ 632 628 624

620-1-'-:---1---1-•-,—

0,0 0,2 0,4 0,6 0,8

Рис. 5. Субоптимальный температурный режим, вычисленный с помощью метода

вариаций в пространстве управления

Решение задачи с целевым функционалом (14) с начальным приближением управления, равным 632 К, совпадает с решением, найденным на основе алгоритма дифференциальной эволюции (Т = 620 К). Однако решение задачи с начальным приближением, равным 644 К, попадает в локальный экстремум. Наименьшая концентрация углекислого газа равна 0,2 моль/л, а оптимальная температура равна 644 К на протяжении всего времени протекания реакции.

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

5. Заключение

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

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

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

6. Благодарности

Исследование выполнено в рамках государственного задания Министерства науки и

высшего образования Российской Федерации (код научной темы FZWU-2023-0002).

7. Литература

[1] Бояринов А.И., Кафаров В.В. Методы оптимизации в химической технологии. М.: Химия, 1975.

[2] Santos LR., Villas-Boas F., Oliveira A.R.L., et al. Optimized choice of parameters in interiorpoint methods for linear programming // Computational Optimization and Applications. 2019. Vol. 73. P. 535-574.

[3] Ziyatdinov N.N., Emel'yanov I.I., Lapteva T.V., et al. Method of Automated Synthesis of Optimal Heat Exchange Network (HEN) Based on the Principle of Fixation of Variables // Theoretical Foundations of Chemical Engineering. 2020. Vol. 54, No. 2. P. 144-162.

[4] Biegler L.T. Integrated Optimization Strategies for Dynamic Process Operations // Theoretical Foundations of Chemical Engineering. 2017. Vol. 51, No. 6. P. 910-927.

[5] Antipina E.V., Antipin A.F., Mustafina S.A. Search for the Optimal Regime Parameters of a Catalytic Process Based on Evolutionary Computations // Theoretical Foundations of Chemical Engineering. 2022. Vol. 56, No. 2. P. 162-169.

[6] Dadebo S.A., Mcauley K.B. Dynamic Optimization of Constrained Chemical Engineering Problems Using Dynamic Programming // Computers & Chemical Engineering. 1995. Vol. 19, Issue 5. P. 513-525.

[7] Островский Г.М., Волин Ю.М. Методы оптимизации сложных химико-технологических схем. М.: Химия, 1970.

[8] Mohamed A.W., Mohamed A.K. Adaptive guided differential evolution algorithm with novel mutation for numerical optimization // International Journal of Machine Learning and Cybernetics. 2019. Vol. 10. P. 253-277.

[9] Xue B., Zhang M., Browne W.N., et al. A Survey on Evolutionary Computation Approaches to Feature Selection // IEEE Transactions on Evolutionary Computation. 2016. Vol. 20, No. 4. P. 606-626.

[10] Рутковская Д., Пилиньский М., Рутковский И. Нейронные сети, генетические алгоритмы и нечеткие системы. М.: Горячая линия-Телеком, 2013.

[11] Карпенко А.П. Эволюционные операторы популяционных алгоритмов глобальной оптимизации // Математика и математическое моделирование. 2018. № 1. C. 59-89.

[12] Антипина Е.В., Мустафина С.А., Антипин А.Ф. Численный алгоритм идентификации кинетической модели химической реакции // Вестник Технологического университета. 2019. Т. 22, № 9. С. 13-17.

[13] Storn R., Price K. Differential Evolution - A Simple and Efficient Heuristic for global Optimization over Continuous Spaces // Journal of Global Optimization. 1997. Vol. 11. P. 341-359.

[14] Антипина Е.В., Мустафина С.А., Антипин А.Ф. Программное обеспечение для автоматизации процесса поиска кинетических параметров химических реакций // Программные продукты и системы. 2020. № 1. С. 125-131.

[15] Григорьев И.В., Михайлова Т.А., Мустафина С.А. О численном алгоритме метода вариаций в пространстве управлений // Фундаментальные исследования. 2015. № 5-2. С. 279-283.

Numerical algorithm for searching the optimal temperature regime of a chemical process on the basis of evolutionary

calculations

Antipina E.V.*, Mustafina S.A.**, Antipin A.F.***

Ufa University of Science and Technology

*

stepashinaev@ya.ru

**

mustafina_sa@mail.ru

***

andrejantipin@ya.ru

Abstract. An algorithm for calculating the optimal temperature profile of a chemical process is developed. The formulation of the problem of optimal control of the chemical process is given. The temperature of the reaction mixture is considered as a control, the values of which are subject to constraints. A finite-dimensional approximating problem is formulated, for the solution of which an algorithm based on the method of differential evolution is constructed. The work of the algorithm is tested on an industrially significant process of phthalic anhydride production. The optimal temperature regimes at which the highest degree of transformation of the starting substance, the highest concentration of the target product and the lowest content of the reaction by-product are achieved are calculated.

Keywords: optimal control problem, optimal temperature regime, chemical process, differential evolution method.

Acknowledgements: This research was funded by the Ministry of Science and Higher Education of

the Russian Federation (scientific code FZWU-2023-0002).

*