CC BY
36
УДК 544.4:004.4 Дата подачи статьи: 21.09.19
Б01: 10.15827/0236-235Х.129.125-131 2020. Т. 33. № 1. С. 125-131
Программное обеспечение для автоматизации процесса поиска кинетических параметров химических реакций
Е.В. Антипина 1, к.ф.-м.н, младший научный сотрудник, е.и.апИрта@р1гЬзи.ги С.А. Мустафина 1, д.ф.-м.н., профессор, тиз1аАпа_эа@таИ-ги А.Ф. Антипин 1, к.т.н, доцент, апйте]апИрт@уа.ги
1 Стерлитамакский филиал Башкирского государственного университета, г. Стерлитамак, Республика Башкортостан, 453103, Россия
В статье описано разработанное авторами ПО для решения задачи поиска кинетических параметров химических реакций.
В условиях современного промышленного производства широко применяются методы математического моделирования химико-технологических процессов, позволяющие на этапе вычислительного эксперимента решать задачи прогнозирования и оптимизации производства. Для решения этих задач необходимо математическое описание химического процесса, которое опирается на построение и исследование кинетической модели химической реакции. Актуальность разработки ПО для решения задачи поиска кинетических параметров химических реакций обусловлена большими вычислительными и временными затратами.
Разработанное ПО применяет метод дифференциальной эволюции и позволяет рассчитать по экспериментальным значениям концентраций веществ значения предэкспоненциальных множителей и энергий активации введенной пользователем реакции. На основе найденных кинетических параметров строится кинетическая модель реакции и решается прямая кинетическая задача. Результаты расчетов выводятся в численном и графическом видах. Для оценки близости полученного решения прямой задачи с опытными данными на одном графике приводятся экспериментальные и численные значения концентраций веществ, а также значения относительных погрешностей концентраций.
Апробация программного комплекса по расчету кинетических параметров проведена на примере реакции получения фталевого ангидрида, в результате сделаны выводы об удовлетворительном согласовании экспериментальных данных со значениями концентраций, вычисленными на основе рассчитанных кинетических параметров.
Программное средство может быть использовано в вычислительных экспериментах научными и производственными лабораториями на этапах проектирования и модернизации химических процессов.
Ключевые слова: программное обеспечение, обратная кинетическая задача, метод дифференциальной эволюции, схема реакции.
При построении математического описания химико-технологических процессов требуется определение динамики взаимодействующих веществ химических реакций. Одним из этапов математической обработки экспериментальных данных является построение и исследование кинетической модели химической реакции, которая представляет собой систему обыкновенных дифференциальных уравнений. Обратная кинетическая задача заключается в определении вида кинетической модели реакции и ее параметров на основе экспериментальных данных [1]. В процессе поиска кинетических параметров химической реакции необходимо минимизировать функционал отклонения экспериментальных данных от расчетных значений, то есть решить оптимизационную задачу.
Однако большинство методов оптимизации, применяемых для решения обратной за-
дачи, зависит от начального приближения поиска решения [2-4]. Одним из методов, позволяющих отыскать глобальный оптимум независимо от начальной точки поиска, является метод дифференциальной эволюции. Это прямой метод оптимизации, так как в ходе его работы не требуется нахождение производных целевой функции, а необходимо определять значения только целевой функции [5, 6].
Расчет кинетических параметров химической реакции является трудоемким процессом, поскольку для вычисления значений целевой функции необходимо неоднократно решать систему дифференциальных уравнений, то есть прямую кинетическую задачу. В связи с этим возникает необходимость в создании программного средства, позволяющего автоматизировать все этапы решения обратной кинетической задачи, что существенно сократит вы-
числительные затраты исследователя при составлении математического описания химических процессов.
Алгоритм поиска кинетических параметров реакции
Пусть кинетическая модель химической реакции описывается системой дифференциальных уравнений:
dx.
-f = f (k, x), i = 1, dt
(1)
где х = (х1, х2, ..., хп) - вектор концентраций веществ, вид функции/определяется схемой реакции; к = (к1, к2, ..., кт) - вектор констант скоростей стадий, причем константа скорости ]-й стадии определяется из уравнения Аррениуса: С Е.
(2)
ki (T ) = k0 , exP
_j_
RT
j = m
где к0]- - предэкспоненциальный множитель; Е] - значение энергии активации]-й стадии; Т-температура; Я - универсальная газовая постоянная [7].
Неизвестными параметрами кинетической модели реакции являются вектор значений предэкспоненциальных множителей к0 = (к01, ко2, ..., кот) и вектор значений энергий активации Е = (Е1, Е2, ..., Ет) стадий реакции.
Требуется найти такой набор значений кинетических параметров ко и Е, при которых функционал отклонения между расчетными и экспериментальными данными принимает минимальное значение:
ß = ZZ| x; - x
(3)
где х^ - значения концентраций веществ, полученных в результате решения системы дифференциальных уравнений (1); хЭ - экспериментальные значения концентраций веществ; к - количество точек эксперимента; п - количество веществ.
Сформулируем алгоритм поиска кинетических констант с помощью метода дифференциальной эволюции. Данный метод имитирует этапы эволюции живых организмов [8, 9]. Аналогом поколений живых существ выступает набор векторов с фиксированным числом координат к] = (к]1, ] ..., к]т),] = 1, ..., соыМРор (соыМРор задает количество векторов в популяции), которые представляют собой кинетические константы скоростей реакции. Представим алгоритм в виде последовательности шагов.
Шаг 1. Создание начальной популяции.
Задаются параметры алгоритма: размер популяции countPop, максимальное количество популяций max Pop, весовой масштабирующий множитель F е [0.4, 1], параметр оператора скрещивания Cross е [0, 1], счетчик итераций i = 0.
Случайным образом генерируется начальная популяция констант скоростей стадий к1, кг, ..., kcount Pop, где kj = (j к/2, ..., km),j = 1, ..., countPop, в пределах пространства поиска.
Устанавливается номер вектора-мишени mish = 1, и вектору-мишени kmish присваивается значение первой константы скорости к1.
Для каждого набора кинетических констант вычисляется значение целевой функции (3) путем решения системы дифференциальных уравнений (1).
Шаг 2. Мутация.
Из текущей популяции случайным образом выбираются два вектора kq, kp так, что q Ф mish, p Ф mish, и третий вектор ki (l Ф mish) с наименьшим значением целевой функции (3). Генерируется новый вектор-мутант по формуле kmut =
= ki + F(kq - kp).
Шаг 3. Кроссовер.
Создается пробный вектор kprob следующим образом. Для каждой координаты вектора kj = = (k/1, k/2, ., kjm), j = 1, count Pop, на интервале [0, 1] генерируется случайное число r (i = 1, ..., m). Если сгенерированное число меньше значения параметра оператора скрещивания (ri < Cross), то kprob i := kmut i, иначе kprob i := := kmish i. Для пробного вектора kprob и вектора-мишени kmish вычисляются значения целевой функции (3).
Шаг 4. Создание новой популяции.
Если Q(kprob) < Q(kmish), то в новую популяцию помещается пробный вектор, иначе - вектор-мишень.
Осуществляется переход к следующей итерации работы алгоритма. Если mish < count_ Pop, то mish := mish + 1 и выполняется переход на шаг 2, иначе - на шаг 5.
Шаг 5. Проверка условия окончания поиска.
Если i < max_Pop, то i := i + 1 и осуществляется переход к шагу 2, иначе - к шагу 6.
Шаг 6. Формирование решения задачи поиска кинетических констант.
Из последней популяции k\, k2, ..., kcount Pop выбирается вектор констант, который соответствует наименьшему значению целевой функции. Этот вектор будет решением задачи поиска кинетических констант стадий реакции.
n
Э
^ mm
На основе вычисленных значений констант скоростей стадий реакции определяются вектор значений предэкспоненциальных множителей ¿0 = (&01, ¿02, ..., ¿0т) и вектор значений энергий активации Е = (Е\, £2, ..., Ет) стадий реакции.
Из уравнения Аррениуса (2) рассчитываются значения энергий активации:
1п к = 1п к„ - .
ЯТ
Если известны значения кинетических констант скоростей стадий реакции при различных значениях температуры к(Т\), к(Т2), ..., кг(Тмт), г = 1, ..., т, то получаем систему уравнений:
ln к. (1 ) = ln kn.---,
'( 1) 0l RT
ln ki (T2 ) = ln koi --
E.
1
RT
(4)
ln k (tnt ) = ln k0i -
E _l_
RT
E.
Обозначим ln k„. = a,--'- = b
R
1
1пк,(Т) = у, -= х], ] = 1, ...,МТ.
}
Тогда система уравнений (4) примет вид
у1 = а + Ъ%1, у2 = а + Ъх2,
уш = а + Ъхш .
Неизвестными параметрами данной системы уравнений являются а и Ь, для вычисления которых методом наименьших квадратов требуется минимизировать функционал:
N1 2
R(a, b) = £((a + bxj )- y, )
(5)
j=1
Значения энергий активаций E и константы ко рассчитываются по формулам:
ко = ea, E = - bR. (6)
Сформулированный алгоритм программно реализован на языке программирования Delphi. Разработанное приложение позволяет по набору экспериментальных данных вычислять кинетические параметры химических реакций.
Описание ПО
Работа программы состоит из нескольких этапов.
Этап 1. Считывание из файла схемы реакции.
С помощью пункта меню «Файл» из текстового файла считывается схема реакции в специальном формате, который включает в себя обозначения реагентов и продуктов реакции. Например, реакция вида X + Х2 ^ Хз должна быть представлена в файле в виде строки: Х1 + Х2 ^ Х3.
Синтаксический анализатор программы на основе схемы реакции выстраивает кинетическую модель, то есть систему обыкновенных дифференциальных уравнений.
Этап 2. Считывание/ввод экспериментальных данных.
Для расчета значения функционала отклонения между расчетными и экспериментальными данными необходимо ввести/считать из файла значения концентраций веществ, полученные опытным путем. Это можно сделать на вкладке «Экспериментальные данные» (рис. 1) либо непосредственно введя в таблицу значения концентраций реагентов, либо открыв данные из текстового файла («Считать из файла»), причем первый столбец значений означает текущее время, остальные - текущие значения концентраций веществ.
Рис. 1. Окно программы на этапе установки параметров оптимизационного алгоритма
Fig. 1. The program window at the stage of setting the parameters of the optimization algorithm
Этап 3. Ввод параметров процесса. В ходе решения оптимизационной задачи многократно решается построенная система дифференциальных уравнений, поэтому пользователю необходимо ввести на вкладке «Параметры процесса» значения начальных концентраций веществ, температуру и время протекания реакции.
Этап 4. Установка параметров алгоритма.
На вкладке «Параметры алгоритма» необходимо ввести такие параметры, как размер популяции, количество популяций, количество констант скоростей стадий, ограничения на константы скоростей реакции, весовой масштабирующий множитель (число из интервала [0.4, 1]), параметр оператора скрещивания (число из интервала [0.4, 1]).
После установки необходимых параметров нужно нажать на кнопку «Расчет к». В результате вычислений определяются значения кинетических констант реакции. Чтобы рассчитать значения предэкспоненциальных множителей и энергий активации стадий, необходимо перейти на вкладку «Кинетические параметры» и нажать на кнопку «Расчет к0 и Е». На основе найденных кинетических параметров решается прямая кинетическая задача, решение ее представляется в численном и графическом видах (вкладки «Решение прямой задачи», «Графики»). На вкладке «Графики» строится график динамики концентраций веществ, рассчитанных по вычисленным кинетическим параметрам и найденных экспериментальным путем. Пользователь может отобразить концентрации всех веществ или части из них. Кроме этого, для каждого вещества вычисляется относительная погрешность значений концентраций, что дает возможность сделать вывод об удовлетворительном согласовании опытных данных с рассчитанными значениями.
Вычислительный эксперимент
Применяя разработанное программное средство, вычислим кинетические параметры реакции получения фталевого ангидрида.
Фталевый ангидрид - один из важнейших видов сырья для производства лакокрасочных материалов, пластификаторов, лекарственных веществ, красителей, присадок к смазочным маслам, ускорителей вулканизации каучука, добавок к реактивным топливам, инсектицидов.
Реакция получения фталевого ангидрида представляется совокупностью следующих стадий [10]: XI ^ Х2 Хг ^ Х4
Х1 ^ Хз (7)
Х1 ^ Х4 Х2 ^ Хз Хз ^ Х5,
где Х1 - нафталин; Х2 - нафтохинон; Хз - фталевый ангидрид; Х4 — углекислый газ; Х5 - ма-леиновый ангидрид.
Кинетическая модель данной реакции представляет собой систему дифференциальных уравнений: (Хх
--Хл ^ Хл к л Хл ,
dt
(8)
dx2
— к^Х^ к2Х2 к^Х^ ,
dt
(Хъ
— к-> Хл + к с х~ к /- Х-,, », 31 52 63'
dt
— к^ Х2 к^ Х^ ,
dt (Х.
- - к,х~ ,
6 3'
где х - концентрация /-го вещества (/ = 1, ..., 5) (моль/л); к - константа скорости'-й стадии реакции (7 = 1, ..., 6) (1/ч).
На основании экспериментальных данных, представленных в [10], для реакции получения фталевого ангидрида рассчитаны значения кинетических констант и энергии активации (рис. 2).
В ходе работы алгоритма система дифференциальных уравнений (8) с начальными условиями решалась методом прогноза и коррекции. Вычисления проводились при следующих параметрах оптимизационного алгоритма: размер популяции - 60, количество популяций - 2 000, количество констант скоростей стадий - 6, весовой масштабирующий множитель - 0.6, параметр оператора скрещивания -0.7. Начальные концентрации веществ задавались следующими значениями (моль/л): Х1(0) = 1, х(0) = 0, / = 2, ..., 5.
Время протекания реакции - 0.6 ч.
На рисунке 2 приведены численные значения кинетических констант и энергий активации стадий реакции, а также представлены графики динамики концентраций веществ Х\ и Х2, вычисленных на основе найденных кинетических параметров и полученных экспериментальным путем при температуре 620 °К.
Заключение
Чтобы определить, согласуются ли экспериментальные данные с рассчитанными значениями концентраций, в программе для каждого вещества рассматриваемой реакции (7) вычислены относительные погрешности нахождения
Рис. 2. Результаты решения обратной кинетической задачи Fig. 2. Results of solving the inverse kinetic problem
значений концентраций веществ: 8(Х1) = = 11,5 %, 8(Х2) = 8,5 %, 8(Х4) = 9 %, 8Х5) = = 11,2 %.
Поскольку отклонение расчетных значений концентраций, полученных на основе вычисленных кинетических параметров, от опытных данных не превышает 12 %, можно сделать вывод об удовлетворительном согласовании с результатом эксперимента.
Таким образом, разработанная программа позволяет рассчитать значения констант скоростей и энергий активации ее стадий. Кинетиче-
ская модель реакции, построенная на основе вычисленных кинетических параметров, достаточно близко описывает динамику концентраций веществ, полученных экспериментальным путем.
Программное средство может применяться в практической деятельности научно-исследовательских лабораторий промышленных предприятий при проведении вычислительных экспериментов разрабатываемых процессов, а также при разработке их математического обеспечения.
Исследование выполнено при финансовой поддержке РФФИ и Правительства Республики Башкортостан в рамках научного проекта № 17-47-020068.
Литература
1. Губайдуллин И.М., Рябов В.В., Тихонова М.В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислительные методы и программирование. 2011. T. 12. № 1. С. 127-135.
2. Исмагилова А.С., Пыжьянова Л.Р., Спивак С.И. Автоматизация метода декомпозиции для анализа информативности кинетических измерений при решении обратных задач // Вестн. Башкирского ун-та. 2016. Т. 21. № 1. С. 9-13.
3. Седова Н.А., Седов В.А. Логико-лингвистическая модель оценки уровня аварийных ситуаций // Современная наука: актуальные проблемы теории и практики. 2016. № 2. С. 65-69.
4. Степашина Е.В. Оптимизация финансовых показателей предприятия на основе нейросетевой модели // Информационные системы и технологии. 2014. № 5. С. 34-42.
5. Canayaz M., Karci A. Cricket behavior based evolutionary computation technique in solving engineering optimization problems. Applied Intelligence, 2016, vol. 44, pp. 362-376.
6. Zaheer H., Pant M., Kumar S., Monakhov O., Monakhova E., Deep K. A new guiding force strategy for differential evolution. Intern. J. of System Assurance Engineering and Management, 2015, vol. 6, pp. 1-14. DOI: 10.1007/s13198-014-0322-6.
7. Антипина Е.В., Антипин А.Ф. Алгоритм расчета оптимальных начальных концентраций веществ химических реакций // Вестн. Технологического ун-та. 2017. Т. 20. № 13. С. 84-87.
8. Karci A. Differential evolution algorithm and its variants. Anatolian J. of Computer Sciences, 2017, vol. 2, no. 1, pp. 10-14.
9. Седова Н.А. Нечеткая продукционная модель первичной оценки опасности столкновения судов // Мир транспорта. 2015. Т. 13. № 2. С. 200-206.
10. Вайтиев В.А., Мустафина С.А. Построение двусторонних оценок решения прямой задачи химической кинетики // Журнал Средневолжского математического об-ва. 2012. Т. 14. № 4. С. 18-25.
Software & Systems Received 21.09.19
DOI: 10.15827/0236-235X.129.125-131 2020, vol. 33, no. 1, pp. 125-131
Software for automation of the search process for kinetic parameters
of chemical reactions
E.V. Antipina 1, Ph.D. (Physics and Mathematics), Junior Researcher, [email protected] S.A. Mustafina 1, Dr.Sc (Physics and Mathematics), Professor, [email protected] A.F. Antipin 1, Ph.D. (Engineering), Associate Professor, [email protected]
1 Sterlitamak Branch of Bashkir State University, Sterlitamak, 453103, Russian Federation
Abstract. The paper describes the software, which the authors developed. It is for solving the problem of finding the kinetic parameters of chemical reactions. In the conditions of modern industrial production, methods of mathematical modeling of chemical-technological processes are widely used, which allow solving the problems of forecasting and optimizing production at the stage of a computational experiment. To solve these problems, we need the mathematical description of the chemical process, which relies on the construction and observation of the kinetic model of the chemical reaction.
The relevance of software development for solving the problem of searching for kinetic parameters of chemical reactions is due to large computational and time costs. The software developed by the authors uses the differential evolution method and allows calculating the values of preexponential factors and activation energies of the user-entered reaction from the experimental values of the concentrations of substances. Based on the found kinetic parameters, the authors construct a kinetic reaction model and solve the direct kinetic problem. The calculation results are in numerical and graphical form. For appraisal the proximity of the obtained solution of the direct problem with the experimental results, there are the experimental and numerical values of the concentrations of substances, as well as the values of the relative errors of the concentrations on one graph.
Software package approbation of the kinetic parameters calculation was the case of the phthalic anhydride reaction. In view of the results, there were implications about the satisfactory agreement of the experimental data with the concentrations calculated because of the calculated kinetic parameters.
We can use the software tool during the computational experiments in scientific and production laboratories at the design and modernization stages of chemical processes.
Keywords: software, inverse kinetic problem, differential evolution method, reaction scheme.
Acknowledgements. The reported study was funded by RFBR according to the research project no. 17-47-020068.
References
1. Gubaydullin I.M., Ryabov V.V., Tikhonova M.V. Application of the global optimization index method to solving inverse problems of chemical kinetics. Numerical Methods and Programming. 2011, vol. 12, no. 1, pp. 127-135 (in Russ.).
2. Ismagilova A.S., Pyzhyanova L.R., Spivak S.I. Automation of decomposition method for the analysis of informativity of kinetic measurements in the inverse problem. Bull. of Bashkir University. 2016, vol. 21, no. 1, pp. 9-13 (in Russ.).
3. Sedova N.A., Sedov V.A. The logical-linguistic model for assessing of the emergency level. Modern Science: Actual Problems of Theory and Practice. 2016, no. 2, pp. 65-69 (in Russ.).
4. Stepashina E.V. Optimization the financial performance of enterprise on neural network model. Information Systems and Technologies. 2014, no. 5, pp. 34-42 (in Russ.).
5. Canayaz M., Karci A. Cricket behavior based evolutionary computation technique in solving engineering optimization problems. Applied Intelligence. 2016, vol. 44, pp. 362-376.
6. Zaheer H., Pant M., Kumar S., Monakhov O., Monakhova E., Deep K. A new guiding force strategy for differential evolution. Int. J. of Syst. Assurance Eng. and Management. 2015, vol. 6, pp. 1-14. DOI: 10.1007/s13198-014-0322-6.
7. Antipina E.V., Antipin A.F. Algorithm for calculating the optimal initial concentrations of substances of chemical reactions. Bulletin of Technological University. 2017, vol. 20, no. 13, pp. 84-87 (in Russ.).
8. Karci A. Differential evolution algorithm and its variants. Anatolian J. of Computer Sciences. 2017, vol. 2, no. 1, pp. 10-14.
9. Sedova N.A. Fuzzy production model for initial evaluating of the risk of collisions. World of Transport and Transportation J. 2015, vol. 13, no. 2, pp. 200-206 (in Russ.).
10. Vaytiev V.A., Mustafina S.A. Creation of bilateral estimates of the solution of a direct problem of chemical kinetics. Middle Volga Mathematical Society J. 2012, vol. 14, no. 4, pp. 18-25 (in Russ.).
Для цитирования
Антипина Е.В., Мустафина С.А., Антипин А.Ф. Программное обеспечение для автоматизации процесса поиска кинетических параметров химических реакций // Программные продукты и системы. 2020. Т. 33. № 1. С. 125-131. DOI: 10.15827/0236-235X.129.125-131.
For citation
Antipina E.V., Mustafina S.A., Antipin A.F. Software for automation of the search process for kinetic parameters of chemical reactions. Software & Systems. 2020, vol. 33, no. 1, pp. 125-131 (in Russ.). DOI: 10.15827/0236-235X. 129.125-131.
