Кинетика и катализ, 2020, T. 61, № 6, стр. 783-788

ВВЕДЕНИЕ

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

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

ТЕОРЕТИЧЕСКАЯ ЧАСТЬ

Рассмотрим произвольную химическую реакцию, протекающую через стадии вида

где A_k – реагенты, k = 1, …, K – номер реагента, a_±ij ≥ 0 − стехиометрические коэффициенты реагента A_k в стадии i = 1, …, s. Динамика такой реакции в открытом безградиентном изотермическом реакторе в рамках закона действующих масс описывается системой K обыкновенных дифференциальных уравнений (ОДУ) где ${{r}_{{ \pm i}}} = {{k}_{{ \pm i}}}{{\Pi }_{k}}A_{k}^{{a \pm ik}}$ − скорости стадий в прямом и обратном направлениях (1/с), k_±i − константы скоростей стадий, q⁰, q и $A_{k}^{0}$ и A_k − соответственно, скорости подачи (1/с) и безразмерные концентрации реагентов на входе и выходе из реактора, A_k(t₀) = $A_{k}^{0}$ − начальные условия (н. у.). Координаты стационарных состояний (с. с.) $A_{k}^{\infty } = A_{k}^{\infty }({{k}_{{ \pm i}}},{{q}^{0}},A_{k}^{0})$ определяются из условий

Если для реакции выполняются линейные стехиометрические законы сохранения (ЛСЗС) вида

где α_jk и C_j – константы, то размерность системы (1) можно уменьшить на число этих ЛСЗС. Как показано в работе [26], точное число ЛСЗС равно где R_s − ранг матрицы A = (a_ik − a_−ik), R_k − ранг по комплексам констант скоростей стадий. Выразим c помощью ЛСЗС число зависимых реагентов N_s через остальные (независимые) и исключим их из системы (1)−(2). Если N_s = 0, то система (1)−(2) не меняется. С учетом этого далее будем считать, что уравнения (1)−(2) включают только независимые реагенты и их число равно К.

Пусть наблюдаемые концентрации всех реагентов характеризуются монотонным изменением во времени. Составим таблицу их значений A_kn от времени t_n с шагом h, где n = 1, 2, 3, … − номер измерения. Выберем из этой таблицы набор из N точек (таких наборов может быть несколько, и можно выбрать любой). Каждый такой набор соответствует некоторому кусочно-линейному приближению (линейному сплайну), заменяющему кривые концентраций вписанными непрерывными ломанными линиями, состоящими из N − 1 отрезков прямых (фрагментов). Выберем опорные (используемые в дальнейших расчетах) значения концентраций реагентов для каждого фрагмента как средние арифметические между их границами

где A_kn и A_{kn − 1} – концентрации реагентов на правой и левой границах n-го фрагмента (индексы 0 и N соответствуют н. у. A_k0 ≡ $A_{{k0}}^{0}$ и c. c. A_kN ≡ A_k∞). Найдем средние скорости изменений концентраций реагентов для каждого фрагмента (можно использовать экспериментальные данные, если они есть) где H_n = N_nh − длина n-го фрагмента, N_n − число интервалов измерений на n-ом фрагменте. Учитывая, что разные фрагменты экспериментальных данных характеризуются разной информативностью, разделим их на три участка: 1) быстрая релаксация (вдали от с. с.) − короткий и очень информативный; 2) обычная релаксация (от конца первого участка до окрестности с. с.) – средний по длине и информативности; 3) медленная релаксация (вблизи с. с., квазистационарность) − длинный и слабо информативный участок. Определим визуально границы (необходимы минимум три точки) наиболее информативных фрагментов и опишем их линейными сплайнами длиной H_n. Подставим выражения (5) и (6) в уравнение (1) и получим систему K × N линейных уравнений для определения 2s неизвестных констант скоростей стадий k_±i реакции (I): где $r_{{ \pm in}}^{*} = {{r}_{{ \pm in}}}(A_{{kn}}^{*})$ – скорости стадий, вычисленные в точке (5); k = 1, 2, …, K, n = 1, 2, …, N. Эта система разрешима только тогда, когда

При K × N = 2s система (8) имеет единственное решение

где Δ ≠ 0 и Δ_±i – главный и вспомогательный определители системы (7) соответственно. Отметим, что вырожденные решения, обусловленные особой структурой кинетических уравнений (см., например, [9, 11]), здесь не рассматриваются. В этом случае все константы скоростей стадий определяются однозначно, но некоторые (или все) из них могут оказаться нефизичными (отрицательными). Условия физичности констант можно записать в виде

При K × N < 2s существует бесконечно много решений, и константы скоростей стадий определяются неоднозначно. Условия физичности аналогичны (10). Иначе, т.е. если условия (8)−(10) не выполняются, система (7) не имеет физичных решений при данном выборе точек для расчетов. В этих случаях необходимо выбрать другой набор экспериментальных точек. Полный перебор всех возможных комбинаций экспериментальных точек позволяет установить интервалы возможных изменений констант скоростей стадий. Если ни один из наборов экспериментальных точек не дает физичных значений констант скоростей стадий, то обратная задача не разрешима и необходимы новые экспериментальные данные или альтернативная стадийная схема реакции.

Пример 1. Пусть реакция А = С + D протекает по двухстадийной схеме

Для нее уравнения (2) запишутся как

Система (1.1)−(1.2) характеризуется двумя ЛСЗС A + B + C = 1 и A + B + D = 1, т.е. в реакции (II) участвуют два зависимых и два независимых реагента. Выберем в качестве независимых реагентов, например, А и С. Исключим с помощью этих ЛСЗС зависимые реагенты, тогда система (1.1)−(1.2) примет вид

Зададим, например, $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 1, $k_{{ + 2}}^{*}$ = 1, $k_{{ - 2}}^{*}$ = 1 (“истинные” значения констант скоростей стадий) и н. у. реакции q = q⁰ = 1, A⁰ = 1, B⁰ = 0, C⁰ = 0, D⁰ = 0. Численно проинтегрируем систему (1.3)−(1.4) с шагом h на интервале t ∈ [0,t_max]. Получим N точек по каждому реагенту и примем их за исходные данные. Выделим границы фрагментов быстрой, средней и медленной релаксаций. Выберем на двух информативных фрагментах две опорные точки (A₁, С₁) и (A₂, С₂), удовлетворяющие вместе с н. у. условиям физичности (10). Подставим эти значения в (1.3)−(1.4) и получим систему линейных по константам уравнений

Результаты решения этой системы без учета шума (ошибки эксперимента) для разного числа измерений и границ информативных фрагментов приведены в табл. 1.

Из табл. 1 следует, что с ростом числа экспериментальных точек N найденные физичные значения констант стабилизируются (метод устойчив). Ошибки метода не превышают погрешностей измерений <5%. С учетом неоднозначности выбора фрагментов сплайнов решением обратной задачи (без учета шума) являются интервалы значений констант, определенных по пяти и более точкам k₊₁ ∈ [0.7438, 0.8975], k₋₁ ∈ [0.4267, 1.0533], k₊₂ ∈ ∈ [1.4612, 1.7183], k₋₂ ∈ [10.3054, 11.1214]. Эти интервалы (кроме последнего) близки к “истинным” значениям констант скоростей стадий.

Для оценки результативности метода варьировали “истинные” значения констант для 7−13 точек. Например, для “истинных” значений $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 2, $k_{{ - 2}}^{*}$ = 0.2 метод дает k₊₁ = 1.0041, k₋₁ = 0.6996, k₊₂ = 3.0096, k₋₂ = 5.5848 (R = 1.21, E = 4.56). Для значений $k_{{ + 1}}^{*}$ = 2, $k_{{ - 1}}^{*}$ = 0.5, $k_{{ + 2}}^{*}$ = 3, $k_{{ - 2}}^{*}$ = 1 получим k₊₁ = 1.7180, k₋₁ = 0.5997, k₊₂ = = 5.6507, k₋₂ = 7.8756 (R = 1.31, E = 4.90). Даже для “жестких” значений $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 10, $k_{{ - 2}}^{*}$ = = 100 получим k₊₁ = 0.9796, k₋₁ = 0.4975, k₊₂ = = 6.0535, k₋₂ = 59.8843 (R = 1.22, E = 4.60). Как видно, метод достаточно результативен. Для оценки помехоустойчивости метода определили влияние уровня шума в исходных данных на погрешности вычислений (табл. 2).

Из табл. 2 следует, что с ростом уровня шума от 1 до 15% найденные значения констант остаются физичными и решениями обратной задачи (с учетом шума) являются интервалы k₊₁ ∈ [0.8975, 1.3238], k₋₁ ∈ [1.0533, 3.3736], k₊₂ ∈ [1.1882, 1.4612], k₋₂ ∈ [1.6837, 10.3054], которые ближе к “истинным” значениям констант скоростей стадий, чем интервалы (без учета шума).

Пример 2. Пусть реакция А = C протекает через три стадии

Для нее уравнения (1) запишутся как

В этой системе уравнений нет ЛСЗС и все реагенты независимы. Зададим “истинные” значения констант $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 1, $k_{{ + 2}}^{*}$ = 1, $k_{{ - 2}}^{*}$ = 1, $k_{{ + 3}}^{*}$ = 1, $k_{{ - 3}}^{*}$ = 1 и н. у. q = q⁰ = 1, A⁰ = 1, B⁰ = 0, C⁰ = 0. Выберем удовлетворяющие условиям (10) две опорные точки (A₁,С₁) и (A₂,С₂) на информативных фрагментах и подставим их значения в уравнения (2.1)–(2.3):

Результаты решения этой системы без учета шума для разного числа измерений и границ информативных фрагментов приведены в табл. 3.

Как видно из этой таблицы, найденные (без учета шума) константы близки к “истинным” значениям и устойчиво находятся по пяти и более точкам в интервалах k₊₁ ∈ [0.8851, 0.9683], k₋₁ ∈ ∈[1.0229, 1.2922], k₊₂ ∈ [0.8972, 0.9433], k₋₂ ∈ ∈ [1.0981, 1.4686], k₊₃ ∈ [1.12555, 1.6052], k₋₃ ∈ [0.8971, 0.9507]. Для оценки помехоустойчивости метода определили влияние уровня шума в исходных данных на погрешности вычислений (табл. 4).

Из табл. 4 видно, что с ростом уровня шума от 1 до 5% найденные константы остаются физичными и решениями обратной задачи (с учетом шума) являются интервалы k₊₁ ∈ [0.3107, 0.9683], k₋₁ ∈ [0.3105, 1.2922], k₊₂ ∈ [0.6752, 0.9433], k₋₂ ∈ ∈ [0.9201, 1.4686], k₊₃ ∈ [1.6052, 3.3809], k₋₃ ∈ [0.9507, 1.5598], которые близки к “истинным” значениям констант скоростей стадий.

Пример 3. В работе [27] было экспериментально установлено, что реакция гидроалюминирования олефинов для алюминийорганического соединения ${\text{HAlBu}}_{2}^{i}$ (диизобутилалюминийгидрид) в присутствии катализатора Cp₂ZrCl₂ протекает по схеме

где A = [Cp₂ZrH₂⋅ ${\text{ClAlBu}}_{2}^{i}$]₂, B = [Cp₂ZrH₂⋅ ${\text{ClAlBu}}_{2}^{i}$], C = ${\text{HAlBu}}_{2}^{i}$, D = [Cp₂ZrH₂· HAlBu₂ⁱ· ${\text{ClAlBu}}_{2}^{i}$]₂. Для этой схемы в работе [27] c помощью оптимизационных методов при −60°С найдены “истинные” значения констант скоростей стадий (пересчитаны в 1/с): k₊₁ = 0.0328, k₋₁ = 0.0438, k₊₂ = 0.0680, k₋₂ = 0.0037, которые позволяют описывать кинетику данной реакции со стандартной ошибкой 3%. Запишем для схемы (IV) уравнения (1):

где r₊₁ = k₊₁A, r₋₁ = k₋₁B², r₊₂ = k₊₂BC, r₋₂ = k₋₂D. Выпишем ЛСЗС 2A + B + D = 2 и C + D = 1, выберем независимые реагенты А и D, выразим через них B и C и исключим их из выражения (3.1):

Проинтегрируем систему (3.2) при истинных значениях констант с шагом h = 30 на интервале t ∈ [0,200] с н. у. A⁰ = 1, B⁰ = D⁰ = 0, C⁰ = 1. Получим N = 7 “экспериментальных” точек, разделим их на два фрагмента, выберем опорные точки в середине этих фрагментов и вычислим концентрации (A₁,D₁), (A₂,D₂) и скорости $A_{1}^{'},D_{1}^{'}$ и $A_{2}^{'},D_{2}^{'}$ в этих точках. Подставим эти значения в выражения (3.2) и получим линейную систему из четырех уравнений для определения четырех констант скоростей стадий:

Результаты решения этой системы приведены в табл. 5.

Как видно из табл. 5, рассчитанные значения констант c учетом уровня шума 20% находятся в интервалах k₊₁ ∈ [0.0157, 0.0221], k₋₁ ∈ [0.0114, 0.0126], k₊₂ ∈ [0.0578, 0.0986], k₋₂ ∈ [0.0018, 0.0054], которые описывают кинетику данной реакции со стандартной ошибкой, не превышающей 3%.

ЗАКЛЮЧЕНИЕ

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