УДК 544.4
Е. В. Антипина, А. Ф. Антипин
АЛГОРИТМ РАСЧЕТА ОПТИМАЛЬНЫХ НАЧАЛЬНЫХ КОНЦЕНТРАЦИЙ ВЕЩЕСТВ
ХИМИЧЕСКИХ РЕАКЦИЙ
Ключевые слова: химическая кинетика; начальные концентрации веществ; планирование эксперимента; метод искусственных иммунных систем.
В статье сформулирован алгоритм поиска оптимального соотношения начальных концентраций веществ на основе метода искусственных иммунных систем. Алгоритм апробирован на промышленно значимом процессе синтеза бензилиденбензиламина, для которого вычислены оптимальные значения начальных концентраций с целью получения максимального выхода продукта реакции.
Keywords: chemical kinetics; initial concentration of substances; experiment planning; a method of the artificial immune .systems.
The article formulated the algorithm finding the optimal ratio of the initial concentrations on the based the method of artificial immune systems. The algorithm was tested for industrially meaningful process of synthesis of benzi-lidenbenzilamin for which was calculated optimum values initial concentrations in order to obtain maximum yield the reaction product.
Введение
Решение большинства задач в химии и химической технологии связано с проведением сложных и дорогостоящих экспериментов. Однако зачастую при исследовании химического процесса целесообразно на этапе вычислительного эксперимента выявить основные его закономерности ввиду сложности состава исходного сырья и невысокого содержания в нем ценных компонентов. Поэтому разработка методов оптимального планирования эксперимента, позволяющих сэкономить время и материальные средства на проведение натурного эксперимента, является актуальной задачей.
Развитие вычислительной техники и информационных технологий, позволяет применять математические методы при планировании эксперимента в химии. Методы математического моделирования позволяют определить оптимальные концентрации реагентов для получения наилучшего выхода конечного продукта и, следовательно, значительно удешевить себестоимость синтеза веществ.
Расчет начальных концентраций реагентов химической реакции производится путем поиска экстремума функции, вычисляющей динамику концентраций веществ в ходе химического опыта. Существенными недостатками большинства численных методов поиска экстремума функций являются трудности в достижении сходимости процесса, которая напрямую зависит от выбора начального приближения.
В настоящее время широко применяются ме-таэвристические методы оптимизации, которые позволяют эффективно находить глобальный оптимум за приемлемое время [1, 2, 3]. Одной из их особенностей является возможность преодоления точки локального экстремума целевой функции в процессе поиска, поэтому по сравнению с классическими эвристиками они позволяют находить более качественные решения.
Среди метаэвристических алгоритмов выделяется метод искусственных иммунных систем, сильной стороной которого является наличие механизма па-
мяти. Механизм памяти позволяет использовать информацию о найденных ранее наилучших решениях, тем самым повышая эффективность глобального поиска. Важной особенностью метода искусственных иммунных систем является независимость решения оптимизационной задачи от начального приближения [4, 5, 6].
Постановка задачи
Сформулируем в общем виде задачу поиска оптимального соотношения начальных концентраций веществ. Математическая модель динамики концентраций веществ во времени есть система обыкновенных дифференциальных уравнений
dx
— = f (t, x(t ),T), dt
(1)
с начальными условиями
X (0) = х0, / = 1,п, (2)
где х(0 = (х1(?), х2^),...,хп(^)) - вектор концентраций реагирующих веществ, t е ^ ] -промежуток времени функционирования системы, Т - температура протекания реакции [4].
В начальный момент времени вещества связаны между собой некоторым соотношением:
xi(0):...: xn(0) = ,
(3)
Требуется найти такое соотношение начальных концентраций веществ (3), при котором обеспечивается экстремум критерия оптимальности п
0(X) = £л,х,^ вхЪ-, (4)
/=1
выражающий максимальный выход целевых продуктов или минимальное содержание примесей (в зависимости от знака коэффициентов Я) в конечный момент времени.
Алгоритм метода искусственных иммунных систем расчета оптимальных начальных концентраций веществ
Метод искусственных иммунных систем имитирует функционирование иммунной системы живого
n •
организма, назначение которой заключается в уничтожении чужеродных тел и совершенствовании борьбы с ними на основе опыта. Рассмотрим основные понятия метода.
Антигеном называется чужеродное вещество, от которого организм пытается защититься посредством антител. Антитело представляет собой вещество, которое распознает и уничтожает антиген. Клеткой памяти называется иммунная клетка, которая накапливает данные о новых антителах, способных распознать антиген, с целью более эффективной работы в случае попадания в организм подобного антигена. Целевая функция оптимизационной задачи называется функцией приспособленности [7, 8].
В качестве функции приспособленности f будем использовать максимальный выход целевого продукта реакции в конечный момент времени протекания реакции. Для этого потребуется решение прямой кинетической задачи при заданном соотношении начальных концентраций исходных веществ, то есть решения системы дифференциальных уравнений (1) с начальными условиями (2).
Модифицируем алгоритм для решения задачи поиска оптимальных начальных концентраций компонентов реагирующей смеси. Для этого в качестве популяции иммунных клеток введем наборы начальных концентраций веществ:
xj (0) = (x J (0), x 2 (0).....x'n (0)), j =
где k - количество иммунных клеток в популяции.
Для расчета значения функции приспособленности Q( x) необходимо знать концентрации веществ, участвующих в реакции, в конечный момент времени. Для этого введем в алгоритм блок решения прямой кинетической задачи.
Задача поиска минимума целевой функции Q(x) сводится к задаче поиска максимума путем замены знака перед функцией на противоположный: Q(x*) = min Q(x) = max [-Q(x)].
Сформулируем данный алгоритм.
Шаг 1. Создание начальной популяции концентраций исходных веществ.
Задать число иммунных клеток в популяции k, количество родительских клеток для селекции r_sel, количество клеток в популяции с наихудшим значением функции приспособленности n_bad, максимальное количество итераций i_max, параметр оператора клонирования n_cl, параметр мутации m, счетчик итераций num=0.
Случайным образом задать значения k клеток начальной популяции xj (0), j = 1, k, причем 0 < xj (0) < 1.
Шаг 2. Вычисление значений функции приспособленности.
Для сформированных наборов xj(0) вычислить значение функции приспособленности Q(xj) путем решения прямой кинетической задачи для каждой иммунной клетки. То есть на данном шаге необходимо численно найти решение системы дифферен-
циальных уравнений (1) с начальными условиями (2) явным или неявным методом. Начальные условия системы (2) представляют собой иммунную клетку.
Шаг 3. Клонирование.
Выбрать из популяции r_sel родительских клеток (множество антител), которым соответствуют наилучшие значения функции приспособленности (4). Для каждой клетки сгенерировать n_cl клонов.
Шаг 4. Мутация.
Для каждой родительской клетки сгенерировать случайные числа l е [0,1], h1 е [0,1 - х/ (0)],
h2 е [0, х/ (0)], i = 1,n. Выполнить мутацию ее клонов по правилу: если l > 0,5, то х? (0) := х? (0) + h1 • m, иначе xf (0) := xf (0) - h2 • m, где s = 1, n _ cl, m - параметр мутации.
Шаг 5. Вычисление значений функции приспособленности.
Для каждого клона-мутанта вычислить значение функции приспособленности (4).
Шаг 6. Селекция и переход к новой популяции.
Для каждой родительской клетки найти клона-мутанта с наименьшим значением функции приспособленности и сравнить ее значение со значением функции приспособленности родительской клетки. В новой популяции оставить родительскую клетку или заменить ее клоном-мутантом с наилучшим значением функции приспособленности.
Шаг 7. Обновление популяции.
В популяции найти n_bad клеток с наименьшим значением функции приспособленности (4). Найденные клетки заменить новыми, случайно сгенерированными на интервале [0,1), и для них рассчитать значение функции приспособленности (4). Увеличить счетчик итераций num=num+1.
Шаг 8. Проверка условий окончания поиска.
Если num<i_max, то перейти к шагу 3, иначе окончить поиск и перейти к шагу 9.
Шаг 9. Выбор оптимальных концентраций веществ из популяции.
В качестве решения задачи поиска оптимальных концентраций веществ выбрать клетку с наибольшим значением функции приспособленности из последней популяции.
Разработанный алгоритм реализован в виде программного средства в среде визуального программирования Delphi.
Вычислительный эксперимент
На основе сформулированного алгоритма рассчитаем оптимальное соотношение начальных концентраций реакции синтеза бензилиденбензиламина с целью получения максимального выхода продукта реакции. Механизм химической реакции синтеза бензилиденбензиламина в присутствии катализатора FeCl3 • 6H2O представляется совокупностью стадий [9]:
(5)
X1 + X2 ^ X3 + X4,
X3 ^ X5 + X6'
X5 + Xi ^ X7 + X8,
X8 + X6 ^ X9' где X! - бензиламин (C7H9N), X2 - четыреххлори-стый углерод (CCl4), X3 - хлорбензиламин (C7H8NCI), X4 - хлороформ (CHCI3), X5 - 1-фенилметанимин (C7H7N), X6 - хлористый водород (HCl), X7 - бензилиденбензиламин (C14H13N), X8 -аммиак (NH3), X9 - хлористый аммоний (NH4Cl).
Скорость каждой стадии определяется согласно закону действующих масс:
а>1 = ^1X1X2, а>2 = ^2 X3, ®3 = кз X5 Xi, ®4 = ^4 X8 Хб,
где X,- - концентрация /-го вещества ( / = 1,9 ) (мольная доля), kj - константа скорости j-й стадии
реакции (j = 1,4 j, рассчитываемая согласно уравне-
(6)
нию Аррениуса ky (T) = k0у • exp
RT
, где Еу -
энергия активации у-й стадии (кДж/моль), Т - температура протекания реакции (К), R - универсальная газовая постоянная (8,31 Дж/(моль-К)).
Матрица стехиометрических коэффициентов веществ (/у) (/ = 1,4, у = 1,9) приведена в табл. 1.
Таблица 1 - Матрица стехиометрических коэффициентов реакции синтеза бензилиденбензила-мина
a>1 а>2 ®3 (й4
X1 -1 0 -1 0
X2 -1 0 0 0
X3 1 -1 0 0
X4 1 0 0 0
X5 0 1 -1 0
Xa 0 1 0 -1
X7 0 0 1 0
X8 0 0 1 -1
X9 0 0 0 1
£ 0 1 0 -1
Поскольку среди элементов последней строки б есть ненулевые, реакция протекает с изменением реакционного объема [10, 11, 12]. Тогда ее кинетическая модель представляет собой систему дифференциальных уравнений:
¿х, = Fi (х,Т) - X/ • Fп (х,Т) с^ N ,
где
Fi =Y.nWi, i = 1'9'
у=1
где
Fn = £Wy ¿Гу,
у=1 /=1 с начальными условиями:
х,- (0) = x0, i = ITS, N(0) = 1, где N - переменный реакционный
(8) объем,
СО:
W = —-— приведенные скорости стадий реакции C0
( j = 1,4 ) (ч-1), С0 - начальная суммарная концентрация веществ (моль/л).
Функции Fn (x,T), F (x,T) (/ = 1,9) с учетом матрицы стехиометрических коэффициентов имеют вид:
F|( x,T) = -Wi( x,T) - И£( x,T), F2( x,T) = -W|( x,T), Fa( x,T) = Wi( x,T) - W2( x,T), F4 (x,T) = Wi( x,T),
F5 (x, T) = W2( x,T) - W3 (x,T), (9)
F6( x,T) = W2( x,T) - W4( x,T),
F7( x,T) = W3( x,T),
F8( x,T) = W3( x,T) - W4( x,T),
F9( x,T) = W4( x,T),
Fn (x,T) = W2( x,T) - IW4 (x,T).
Исходными веществами реакции являются бензиламин (X1) и четыреххлористый углерод (X2), концентрации которых связаны соотношением: xi(0) + x2(0) = 1.
Для остальных веществ x(- (0) = 0, / = 3,9. Поскольку целевым веществом схемы реакции является бензилиденбензиламин (X7), тогда задача оптимизации формулируется следующим образом: найти такое соотношение начальных концентраций x1(0): x2 (0) исходных веществ X1 и X2, при которых в конечный момент времени протекания реакции достигается максимальный выход бензилиден-бензиламина X7, то есть
Q(x) = x7(t2) ^ max. Решение поставленной задачи получено при температуре T = 23 °C и времени протекания реакции t2 = 4 ч. Решение системы (7) с начальными условиями (8) получено неявным методом Эйлера.
В результате расчета получено оптимальное соотношение начальных концентраций исходных веществ:
x1 (0): x2(0) = 0,51:0,49. (10) При этом максимальный выход продукта реакции (бензилиденбензиламина) составляет 0,013 мольных долей.
На рис. 1 изображена динамика концентрации
бензилиденбензиламина при начальных концентрациях.
оптимальных
Рис. 1 - Динамика оптимальной концентрации бензилиденбензиламина
Результаты и обсуждение
Рассмотрим решение прямой кинетической задачи для схемы реакции (5) при произвольных наборах начальных концентраций веществ. Как видно из табл. 2, максимальное значение концентрации бензилиденбензиламина наблюдается при соотношении исходных веществ (10), что подтверждает корректную работу построенного алгоритма.
Таблица 2 - Зависимость концентрации бензилиденбензиламина Х7 от соотношения исходных
веществ Х1 и Х2
X1, мол. доля X2, мол. доля X7, мол. доля
0,51 0,49 0,013
0,1 0,9 0,004
0,2 0,8 0,008
0,3 0,7 0,010
0,4 0,6 0,012
0,7 0,3 0,011
0,8 0,2 0,009
0,9 0,1 0,005
Таким образом, разработанный алгоритм поиска оптимальных начальных концентраций исходных веществ позволяет на этапе вычислительного эксперимента решить задачу планирования эксперимента в химии. При этом будет найдено решение оптимизационной задачи при любом наборе значений начальных концентраций веществ, поскольку работа метода искусственных иммунных систем не зависит от начального приближения.
Литература
1. Антипина Е.В., Антипин А.Ф. Применение интеллектуальных технологий для анализа многомерных данных // Молодой ученый. 2014. № 19. С. 172-175.
2. Степашина Е.В., Байтимерова А.И., Мустафина С.А. Программный комплекс автоматизации процедуры уточнения механизма химической реакции на основе DRGEP метода // Башкирский химический журнал. 2011. Т. 18. № 3. С. 112-115.
3. Степашина Е.В., Мустафин Т.И. Алгоритм идентификации математической модели редуцированной схемы реакции // Путь науки. 2014. № 8 (8). С. 33-35.
4. Антипин А.Ф. Способ анализа программного кода автоматизированной системы управления технологическими процессами // Автоматизация, телемеханизация и связь в нефтяной промышленности. 2013. № 10. С. 2125.
5. Антипин А.Ф. О проверке программ автоматизированных систем управления // Современная техника и технологии. 2015. № 1. С. 62-66.
6. Мустафина С.А., Степашина Е.В. Редукция кинетических схем сложных химических процессов на основе теоретико-графового подхода // Вестник Казанского технологического университета. 2014. Т. 17. № 10. С. 17-20.
7. Пантелеев А.В., Метлицкая Д.В. Применение метода искусственных иммунных систем в задачах поиска условного экстремума функций // Научный вестник МГТУ ГА. 2012. № 184. С. 54-61.
8. Степашина Е.В. Редуцирование схемы химической реакции // Математические методы в технике и технологиях - ММТТ. 2014. № 8 (67). С. 266-268.
9. Ахметов И.В. Разработка кинетических моделей реакций синтеза ароматических и гетероциклических соединений на основе многоядерных вычислительных систем // Параллельные вычислительные технологии 2013: труды международной научной конференции (г. Челябинск, 01-05 апреля 2013 г.). Челябинск, 2013. С. 268-277.
10. Степашина Е.В., Байтимерова А.И., Мустафина С.А. О свойствах решений задач моделирования каталитических процессов с переменным реакционным объемом // Журнал Средневолжского математического общества. 2010. Т. 12. № 3. С. 122-128.
11. Антипина Е.В., Антипин А.Ф. О существовании решения математической модели химического процесса в реакторе идеального смешения // Фундаментальные и прикладные исследования в современном мире. 2015. № 11-3. С. 32-34.
12. Байтимерова А.И., Степашина Е.В., Мустафина С.А. Математическая модель процесса в РИС на двудольном графе // Обозрение прикладной и промышленной математики. 2010. Т. 17. № 3. С. 462.
© Е. В. Антипина - к.ф.-м.н., доцент каф. математического моделирования Стерлитамакского филиала «Башкирский государственный университет», [email protected]; А. Ф. Антипин - к.т.н., доцент каф. прикладной информатики и программирования Стерлитамакского филиала «Башкирский государственный университет», [email protected].
© E. V. Antipina - Candidate of Physical and Mathematical Sciences, Associate Professor of the Department of Mathematical Modeling, Sterlitamak Branch of Bashkir State University, [email protected]; A. F. Antipin - Candidate of Technical Sciences, Associate Professor of the Department of Applied Informatics and Programming, Sterlitamak Branch of Bashkir State University, [email protected].