Поиск областей неопределенности кинетических параметров математических моделей химической кинетики на основе интервальных вычислений
Автор: Вайтиев Владимир Анатольевич, Мустафина Светлана Анатольевна
Рубрика: Программирование
Статья в выпуске: 2 т.7, 2014 года.
Бесплатный доступ
Проведено исследование влияния неопределенности в кинетических параметрах на результаты решения прямой и обратной задач химической кинетики. Кинетические данные представлены интервалами и рассматриваются как объекты интервального анализа. Разработан алгоритм решения прямой задачи методом интервального анализа чувствительности и алгоритм решения обратной задачи по вычислению областей неопределенности параметров. Проведен вычислительный эксперимент по поиску интервального вектора параметров на примере промышленно значимой реакции. Показано, что интервальное решение прямой задачи, полученное для математической модели реакции, удовлетворяет заданному предельно допустимому значению погрешности в экспериментальных замерах концентраций участвующих веществ.
Интервальный анализ, метод анализа чувствительности, оптимальные границы двустороннего решения, кинетическая модель, константы скорости, прямая и обратная кинетические задачи, интервалы неопределенности кинетических параметров
Короткий адрес: https://sciup.org/147159269
IDR: 147159269 | DOI: 10.14529/mmp140209
Текст научной статьи Поиск областей неопределенности кинетических параметров математических моделей химической кинетики на основе интервальных вычислений
Изучение механизмов сложных химических реакций на. сегодняшний день представляет огромный практический интерес. Кинетическая модель является основой математического моделирования химических реакций. Для построения кинетической модели необходимо знать участвующие в химической реакции вещества, и элементарные стадии. Традиционно при моделировании химических реакций принято использовать средние значения кинетических параметров. Моделирование химических реакций на. основе средних значений параметров не позволяет гарантировать такой режим функционирования, который может возникать в процессе эксплуатации [1]. На самом деле кинетические параметры находятся в некотором интервале возможных значений, так как определены по экспериментальным данным в ходе решения обратной задачи идентификации механизмов сложных химических реакций.
Математическое описание химических реакций представляет собой систему дифференциальных, как правило, нелинейных уравнений с заданными начальными условиями, число неизвестных которой равно числу участвующих в реакции веществ. В то же время непосредственному измерению доступна, только часть из этих веществ. Поэтому возникает обратная задача, определения параметров системы дифференциальных уравнений, воспроизводящих часть ее решений. Следствием недостаточной информативности может стать неединственность решения обратной задачи. Кинетические измерения задаются внутри некоторого интервала. точности, определяемого величиной погрешности измерений, к которым, в том числе, относятся экспериментальные концентрации в определенные моменты времени протекания реакции. Решением обратной задачи определения кинетических параметров становится некоторая область, вариация кинетических констант внутри которой сохраняет требуемое качество описания измерений.
Определить доверительный интервал в задаче математической обработки эксперимента можно методом максимального правдоподобия. В этом случае необходимо знание статистического закона распределения погрешности измерений, которое для реальных систем, как правило, отсутствует. Другая постановка расчета областей неопределенности сделана Л.В. Канторовичем [2]. Она не требует знания информации о статистических свойствах распределения погрешности измерений. Необходимо только знание величины предельно допустимой погрешности эксперимента, информация о которой чаще всего присутствует у экспериментатора. Поиск интервальных значений кинетических параметров, образующих искомую область, предполагает достаточно большой объем вычислений. Для повышения эффективности решения данной задачи наряду с высокопроизводительными параллельными вычислительными системами целесообразно использовать понятия и методы интервального анализа [3–6].
1. Вывод математической модели
В основе составления кинетических уравнений для сложных процессов лежит положение о независимом протекании элементарных реакций [7]. В таких процессах одно и то же вещество X i может принимать участие в качестве исходного или конечного вещества в нескольких элементарных стадиях, т.е. сложная реакция может быть записана в виде:
n
^2 V ij X i = 0, j = 1,m, (1)
i =1
где X i – вещества, участвующие в реакции, ν ij – стехиометрический коэффициент при компоненте X i в j-й стадии, m – число элементарных стадий. Для элементарных реакций значениями ν ij могут быть только 0, 1, 2 и 3 (редко). Для итоговых уравнений таких ограничений не существует. В сложных реакциях стехиометрические коэффициенты могут быть дробные.
Скорость реакции ω i по i-му компоненту выражается изменением количества этого компонента N i в единицу времени в единице реакционного пространства. По отношению к исходным веществам такую характеристику называют скоростью превращения или расходования, по отношению к продуктам – скоростью образования.
В закрытой системе скорость гомогенных реакций –
_ 1 dN i
= Vp dt скорость гетерогенных реакций –
_ 1 dN i ^ i = FPP
где V p – объем реакционного пространства, F – площадь поверхности раздела взаимодействующих фаз.
Скорость гомогенных реакций, выраженная через концентрацию i-го компонента Ci посредством соотношения N = VC, (V — объем реакционной смеси), согласно (2), можно записать в виде:
^ i =
1 VdC i + C i dV
V p dt
Если объем реакционной смеси при протекании химической реакции практически не изменяется (dV ~ 0), то V p = V, и выражение (4) принимает вид:
w i =
V dC i dt .
Следует отметить, что это определение широко распространено в научной литературе, но не может служить в качестве общего определения понятия скорости, так как часто ока-V dC зывается, что величина dt i связана не только с числом актов химического превращения, но и с тем, по какому закону изменяется объем системы. В то же время, это изменение может осуществляться произвольным образом. На практике с протеканием газовых реакций при переменном объеме приходится сталкиваться в тех случаях, когда в реакции изменяется число молекул, а давление в системе поддерживается постоянным. Поэтому при разработке математического описания сложного процесса необходимо учитывать изменение числа молей реакционной смеси (или реакционного объема) в ходе протекания химических реакций [8]. Чтобы определить, существует ли такое изменение в ходе реакций, построим матрицу стехиометрических коэффициентов Г, элементами которой являются стехиометрические коэффициенты участвующих в реакции веществ νij . Ее строки соответствуют элементарным реакциям, а столбцы - веществам. В общем виде стехиометрическую матрицу Г можно представить:
ν 11 ν 12 ... ν 1 m
Г _ V 21 V 22 ... V 2 m
* * * * * * * * * * * *
ν n 1 ν n 2 ... ν nm
Размерность стехиометрической матрицы n х m, где n - число веществ, m - число элементарных реакций. Стехиометрические коэффициенты веществ ν ij , исходных для данной реакции, входят со знаком « — » , продуктов реакции - со знаком « + » . Если вещество не участвует в реакции, его стехиометрический коэффициент принимается равным нулю.
Далее с помощью матрицы Г находят компоненты вектора 5 j по следующему правилу
n
5 j _ ^V ij , j = 1,m. (6)
i =1
Если хотя бы одна из компонент вектора δ j отлична от нуля, то реакция протекает с изменением числа молей.
Согласно закону сохранения массы (2), суммарный материальный баланс для варианта, когда суммарная концентрация C _ ^2П =1 C i изменяется во времени, имеет вид:
dC i d ( Cx i ) q
_ - _ —r' i _1'n (7)
j =1
с начальными условиями: X i (0) _ x 0 , где x i _ CC - концентрация i-го компонента в мольных долях, r j – скорость j-й реакции, определяемая согласно закону действующих масс. Общий вид для элементарной химической реакции:
P r _ knC A l ' l =1
где r – скорость реакции, k – константа скорости реакции, C A α l l – концентрация реагента A l в степени, равной стехиометрическому коэффициенту α l в уравнении реакции, P – количество исходных реагентов.
Замыкает систему уравнений (7) условие нормировки по компонентам реакционной среды:
n
^ x i _ 1. (8)
i =1
Начальная мольная плотность жидкости (начальная суммарная концентрация) C o постоянна при любых температурах. Тогда можно разделить (7) на C o . Получим систему дифференциальных уравнений
^Nx i ) ТТЛ • т
--dt-- = 2^vij Wj, i = 1,n, j=1
где N = С - относительное изменение числа молей реакционной среды, W j = j - приведенные скорости химических реакций. Обозначим
Fi = it Vij Wj, i = 1n(10)
j =i
Преобразуем левую часть системы уравнений (9), применив правило дифференцирования произведения функций и условие нормировки (8). Получим
'. = fi = F' ~ xi Fn, i = ГП, dN = fn+1 = Fn dt Ndt с начальными условиями:
xi(0) = x0, i = 1,n, N (0) = 1.(12)
Функции F n , F i выписываются с учетом матрицы стехиометрических коэффициентов.
Система уравнений (11) с начальными условиями (12) является математической моделью сложной (многостадийной) реакции, учитывающей изменение числа молей в ходе ее проведения.
-
2. Постановка обратной задачи по вычислению областей
неопределенности кинетических параметров
Классическая постановка обратных задач состоит в поиске констант скоростей реакции, минимизирующих критерий соответствия расчета состава многокомпонентной реагирующей смеси экспериментальным данным:
ns
E(k) = E Eix Rj - x Ejh (13)
i =1 j =1
где x i R j – расчетные значения концентраций веществ, участвующих в реакции; x i E j – экспериментальные данные; i – номер компоненты, j – номер замера; k – вектор кинетических констант; n – число веществ, концентрация которых поддается измерению; s – число замеров каждой из компонент.
Решение обратной задачи неоднозначно. Обычно это не число (точка в пространстве констант), а область, возникающая, по крайней мере, по двум причинам: 1) избыточность схемы (большое количество констант); 2) погрешность измерений. Эту область принято называть областью неопределенности. В таком случае целесообразно решать обратную задачу по определению некоторых интервалов, образующих область неопределенности, любая вариация кинетических констант скоростей внутри которых сохраняет качество описания измерений, устанавливаемое величиной, предельно допустимой погрешности эксперимента.
Применим основные понятия и методы теории интервального анализа [3, 4, 5] к решению обратной задачи по определению искомой области. Основная идея интервального анализа состоит в замене арифметических операций и вещественных функций над вещественными числами интервальными операциями и функциями, преобразующими интервалы, содержащие эти числа. Ценность интервальных решений заключается в том, что они содержат точные решения исходных задач.
Пусть γ ij – величина допустимой погрешности измерения концентрации i-го компонента в j-й замер. Тогда результатом решения обратной задачи по вычислению области неопределенности будут некоторые интервальные значения кинетических констант скоростей:
k 1 = [k 1 ; k 1 ], ..., k m = [k m ; k m ], (14)
где k j , k j , j = 1,m есть нижние и верхние границы интервалов соответственно. Суть (14) состоит в том, что при решении прямой задачи (11)-(12) для V k 1 G k i , V k 2 £ k 2 , ..., V k m £ k m могут быть получены расчетные концентрации, удовлетворяющие:
x 1j £ [x 1 j Y 1 j ; x 1j + Y 1 j ] x 1 j , ’”’ x nj £ [x nj Y nj ; x nj + Y nj ] x nj , j 1, s. (15)
Как известно, решение обратной задачи предполагает многократное обращение к решению прямой задачи. Тогда решение системы (11)-(12) с интервальными параметрами (14) приведет к интервальным значениям расчетных концентраций:
x Rj = [x Rj ; x R j 1, -, x R = [x R j ; x Rj 1, j =TTs- (16)
Тогда требуется найти такие интервальные значения кинетических констант скоростей реакции (14), что ns
Q( k ) = ^^ q ( x R , x E ) ^ min, (17)
i=1 j=1
где q(x R , x E ) = max {| x R — x E | , | x R — x E |} - расстояние между двумя интервалами, k -интервальный вектор кинетических констант.
-
3. Алгоритм решения прямой задачи методом интервального анализа чувствительности
Для решения задачи (11) – (12) в условиях (14) будем использовать модифицированный метод интервального анализа чувствительности, ранее описанный в работе [6]. Основная идея метода заключается в анализе частных производных решения по параметрам. Для его реализации используется аппарат интервального анализа. Запишем задачу (11) – (12) в виде системы обыкновенных дифференциальных уравнений:
dx- = f i (t,x,k), i = 1,n, n = q + 1,
где q – количество веществ, вступающих в реакцию, xi(0) = xoi, i = 1,n — 1; xn(0) = xon = N (0) = 1. (19)
Будем искать решение системы (18) с начальными условиями (19) в виде x = ( x (1) , x (2) ,..., x (n) ), где x (i) = [x i ,x i ], x i - нижняя и x i - верхняя границы соотвествующего интервала. Аналогичное представление характерно для вектора констант скоростей (14).
Пусть требуется оценить x (i) - верхнюю границу решения по i-й координате. Для ее оценки система уравнений (18) с начальными условиями (19) примет вид:
′ x = f(t,x,k), k £ k, x(0) = xo £ xo- (20)
Интервальные функции x kj (t), x j (t) можно определить, одновременно решая (18)-(19)
и систему вида: |
x kj' = ^ df i (t,x,k)x kj + f(t,x,k), (23) i =i XXl j x k = 0, i = 1,n, j = 1,m, (24) x 0 j = ^ df i (t, x, k)x ij' , (25) l =1 x l 0 x S 1, при i = j , • • i x j = O ij = < . i,j = 1,n- (26) [0, при i = j; |
Совместное решение (18) – (19), (23) – (26), равно как и системы (20), можно получить приближенным явным методом Рунге–Кутты 4-го порядка точности. При этом решение системы дифференциальных уравнений проводится на основе алгоритма нахождения интервальных расширений для полиномиальных функций.
Оценка нижней границы решения по i-й координате x (i) проводится аналогично с соответствующими корректировками, вытекающими из соображений достижения минимума функции на границах получаемых интервалов.
Метод интервального анализа чувствительности дает возможность определения такого момента времени t opt > 0 [6], до которого найденные границы множества решений (18)-(19) могут считаться оптимальными.
Существуют иные подходы к построению двустороннего решения (18) - (19). Так, в [10] на примере реакции с постоянным реакционным объемом рассмотрен метод, основанный на понятиях монотонности, изотонности функций по переменным и приводящий к представлению исходной системы в виде двух независимых подсистем, параллельное решение одной из которых дает нижнюю, другой – верхнюю границы интервального решения (18) – (19).
4. Алгоритм интервального решения обратной задачи по определению констант скоростей
Процедура решения обратной задачи по вычислению областей неопределенности кинетических констант состоит в поиске интервальных значений констант скоростей стадий, минимизирующих функционал (17). Для минимизации функционала (17) будем использовать модифицированный под интервальные вычисления метод Хука – Дживса, который представляет собой комбинацию исследующего поиска с циклическим изменением переменных и ускоряющего поиска по образцу. Сформулируем алгоритм поиска кинетических констант по шагам.
Шаг 1. Ввод исходных данных: t - время протекания реакции, x (0) — вектор начальных концентраций, x E – матрица интервальных значений экспериментальных концентраций в фиксированные моменты времени, ε – минимальное значение функционала (17) (критерий остановки поиска), h j – шаг вариации по границам каждой координаты (константы скорости), e j – минимально допустимое значение шага вариации по границам каждой координаты в ходе исследующего поиска, imax – максимальное количество итераций поиска, k j down , k j up – величины, образующие область поиска (нижний и верхний предел соответственно); j = 1,m, m - количество констант скоростей.
Шаг 2. Задается стартовая точка (начальное приближение) k j = ( k j , k 2 ,..., k m ), где k j = [k j ; k j ]. В качестве начального приближения может быть выбран вектор вырожденных интервалов. Номер текущей координаты s = 1. Номер итерации поиска i = 0. Вычисляется значение функционала Q( k j ).
Шаг 3. Проводится решение прямой задачи при четырех наборах кинетических параметров:
kj+ = (kj, k2>..., [kj + hs; Й>-> km) при kj + hs < kup, ks— = (kj, k2>-> [kj — hs; Й>-> km) при kj - hs > kdown, kj+ = (kj, k2,..., [kj; k1 + hs],..., km) при kj + hs < ku, ks- = (kj, k2,..., [kj; kj - hs],..., km) при kj - hs > kdown.
Варьирование границ на данном шаге алгоритма также сопровождается контролем выполнения условия ≪ непересечения ≫ верхней и нижней границы текущей координаты. Вычисляются значения функционала Q( k j + ), Q( k j_ ), Q( k j+ ), Q( k j^ ) при заданных наборах параметров.
Шаг 4. Если Q( k j ) < Q( k j + ), Q( k j ) < Q( k j- ), Q( k j ) < Q( k s + ), Q( k j ) < Q( k s — ), то шаг h s уменьшается с переходом на Шаг 3. Уменьшение h s на этом шаге алгоритма ограничивается минимально допустимым значением e s . При отсутствии возможности уменьшать шаг далее текущая координата не изменяется. В противном случае k 1 заменяется на один из векторов k s 1 + , k s 1 - , k s 1 + , k s 1 - , соответствующий минимальному из значений функционалов.
Шаг 5. Если s < m, то s = s + 1 и переход на Шаг 3. Иначе путем исследующего поиска по границам всех координат (Шаг 2 – Шаг 4) получаем новый набор констант k 2 = ( k j , k 2 ,..., k m ). Если условия окончания алгоритма Q( k ) < е или количество итераций поиска, равное установленному максимальному значению i, не выполняются, то переход на Шаг 6.
Шаг 6. Начинается поиск по образцу. Проводится расчет нового набора констант по формуле k 3 = k j + c ( k 2 - k j ), где c – параметр алгоритма. Обычно выбирается равным 2. В случае интервальных вычислений c = [2; 2] - вырожденный интервал.
Шаг 7. Исследующий поиск (Шаг 2 – Шаг 4) для нового набора констант k 3 , за исключением того, что шаги вариации по константам h s на этом этапе метода не уменьшаются. В результате получается набор констант k 4 = ( k j , k 2 ,..., k m ).
Шаг 8. Если вектор k 4 отличен от k 3 , то осуществляем замену k j = k 2 , k 2 = k 4 , i = i+1 и переход на Шаг 6. В противном случае k j = k 2 , i = i + 1 и переход на Шаг 3 к исследующему поиску.
5. Вычислительный эксперимент
Проведем вычислительный эксперимент по решению обратной задачи определения области неопределенности констант скоростей на примере промышленно значимой реакции олигомеризации α-метилстирола. Продукты реакции являются ценным нефтехимическим сырьем и используются в качестве диэлектрических жидкостей, основы смазочных масел, пластификаторов каучуков, реактивного топлива, изоляционного материала, стойкого к радиолизу теплоносителя [11]. Введем обозначения для веществ, участвующих в реакции: X i -α-метилстирол; X 2 – 4-метил-2,4 дифенилпентен-1; X 3 – 4-метил-2,4 дифенилпентен-2; X 4 – 1,1,3-триметил-3-фенилиндан; X 5 – тримеры. Реакция является многостадийной и протекает в условиях изменения количества молей реакционной среды. Экспериментальные данные получены в лаборатории приготовления катализаторов Института нефтехимии и катализа РАН [12].
В результате решения обратной задачи был найден набор параметров k 1 = (k 1 ,..., k]^), при котором наблюдается наименьшее отклонение между расчетными и экспериментальными значениями концентрации веществ, участвующих в реакции. В табл. 1 представлены начальное приближение поиска области неопределенности и интервальный вектор констант скоростей, образующий искомую область.
Таблица 1
Константы скоростей, полученные в ходе решения обратной задачи
Константа скорости |
Начальное приближение k 1 |
Интервальное решение k ⋆ |
k 1 |
1,19575 |
[ 1,00825; 1,38325 ] |
k 2 |
0,12781 |
[ 0,10613; 0,14656 ] |
k 3 |
0,19794 |
[ 0,17294; 0,23525 ] |
k 4 |
0,02630 |
[ 0,02318; 0,02552 ] |
k 5 |
0,06974 |
[ 0,06896; 0,06974 ] |
k 6 |
0,49184 |
[ 0,47934; 0,50590 ] |
k 7 |
0,00001 |
[ 0,00020; 0,00040 ] |
k 8 |
0,41116 |
[ 0,40178; 0,41116 ] |
k 9 |
0,07412 |
[ 0,07295; 0,07373 ] |
k 10 |
0,00001 |
[ 0,00001; 0,00176 ] |
k 11 |
0,00649 |
[ 0,00478; 0,00493 ] |
k 12 |
0,00001 |
[ 0,00001; 0,00001 ] |
Интервальное решение обратной задачи получено в ходе минимизации критерия отклонения расчетных и экспериментальных значений концентраций x R и x E (i = 1, 5 - количество веществ, j = 1, 7 - количество замеров), представленных в табл. 2. Компоненты интервальной матрицы экспериментальных концентраций представляют собой интервалы, содержащие 10%-ную погрешность от значений концентраций, соответствующих решению прямой задачи для набора параметорв k 1 .
На рис. 1 – рис. 3 представлена динамика изменения концентраций трех из пяти веществ, участвующих в реакции: точками • обозначены данные эксперимента, □ - верхние и нижние границы интервальных компонентов матрицы экспериментальных данных x i E j , сплошной линией – результаты решения прямой задачи при наборе констант скоростей k 1 ,
Таблица 2

Рис. 2 . Расчетные и экспериментальные концентрации X 2

Рис. 3 . Расчетные и экспериментальные концентрации X 4
Список литературы Поиск областей неопределенности кинетических параметров математических моделей химической кинетики на основе интервальных вычислений
- Иремадзе, Э.О. Неопределенность в кинетических константах и расчет оптимальной температуры/Э.О. Иремадзе, С.А. Мустафина, С.И. Спивак//Математическое моделирование. -2000. -Т. 12, № 3. -С. 21.
- Канторович, Л.В. О некоторых новых подходах к вычислительным методам и обработке наблюдений/Л.В. Канторович//Сибирский математический журнал. -1962. -Т. 3, № 5. -С. 701-709.
- Калмыков, C.А. Методы интервального анализа/C.А. Калмыков, Ю.И. Шокин, З.Х. Юлдашев. -Новосибирск: Наука, 1986. -222 с.
- Шокин, Ю.И. Интервальный анализ/Ю.И. Шокин. -Новосибирск: Наука, 1981. -112 с.
- Шарый, С.П. Конечномерный интервальный анализ: монография/С.П. Шарый. -М.: Изд-во XYZ, 2013. -726 с.
- Добронец, Б.С. Интервальная математика/Б.С. Добронец. -Красноярск: КГУ, 2004. -219 с.
- Эммануэль, Е.М. Курс химической кинетики/Е.М. Эммануэль, Д.Г. Кнорре. -М.: Высш. шк., 1984. -463 с.
- Мустафина, С.А. Расчет оптимальной температуры в условиях неопределенности по кинетическим константам/С.А. Мустафина, С.И. Спивак//Башкирский химический журнал. -1999. -Т. 6, № 1. -С. 64.
- Добронец, Б.С. Приложения интервального анализа чувствительности/Б.С. Добронец, Е.Л. Рощина//Вычислительные технологии. -2002. -Т. 7, № 1. -С. 75-82.
- Вайтиев, В.А. Численное исследование процессов с постоянным и переменным реакционным объемом в условиях неопределенности кинетических данных/В.А. Вайтиев, С.А. Мустафина//Башкирский химический журнал. -2013. -Т. 20, № 2. -С. 45-48.
- Оптимальные технологические решения для каталитических процессов и реакторов/С.А. Мустафина, Ю.А. Валиева, Р.С. Давлетшин, А.В. Балаев, С.И. Спивак//Кинетика и катализ. -2005. -Т. 46, № 5. -С. 749-756.
- Байтимерова, А.И. Поиск оптимального управления в каскаде реакторов для процессов с переменным реакционным объемом/А.И. Байтимерова, С.А. Мустафина, С.И. Спивак//Системы управления и информационные технологии. -2008. -№ 2(32). -С. 38-42.