Химическая физика, 2021, T. 40, № 6, стр. 47-53

ВВЕДЕНИЕ

Обратная задача химической кинетики, связанная с идентификацией параметров кинетических моделей, относится к некорректным задачам математической физики и химии. Сложность их решения обусловлена погрешностями экспериментальных данных и неоднозначностью решений. Математической базой решения некорректных задач стал метод регуляризации А.Н. Тихонова (1965), который позволяет найти приближенное решение, близкое к точному [1, 2]. Физико-химические исследования обратных задач проводили В.Н. Писаренко и А.Г. Погорелов (1969), М.Г. Слинько (1972), С.И. Спивак (1972), Г.С. Яблонский (1977), В.Н. Лукашенок (1979), В.Г. Горский (1981) и другие. В химической кинетике обратная задача связана с определением констант скоростей элементарных стадий заданного (предполагаемого) механизма по различным экспериментальным данным исследуемой реакции. При ее решении, как правило, используются сложные оптимизационные алгоритмы, которые позволяют определить интервалы изменения значений только некоторых комплексов констант скоростей стадий [3–32 ]. Ниже приведен метод решения обратной задачи, основанный на использовании данных нескольких (двух и более) стационарных экспериментов с разными начальными условиями (мультиэкспериментов) [33–39]. Метод позволяет определять с хорошей точностью точечные физичные значения и интервалы возможного изменения констант скоростей всех элементарных стадий предполагаемого механизма реакции без использования оптимизационных алгоритмов. Показана эффективность применения метода для ряда многостадийных реакций, протекающих в изотермическом безградиентном реакторе, с учетом ошибки определения экспериментальных данных.

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

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

где A_k – реагенты; k = 1, …, K – номер реагента; a_±ik ≥ 0 – стехиометрические коэффициенты реагента A_k в стадии i = 1, …, s. Стационарные режимы такой реакции в изотермическом безградиентном реакторе описываются системой K нелинейных алгебраических уравнений:

где r_±i = ${{k}_{{ \pm i}}}\prod\nolimits_k {A_{k}^{{a \pm ik}}} $ – стационарные скорости стадий в прямом и обратном направлениях (1/с), k_±i – константы скоростей стадий; q⁰, q и $A_{k}^{0},$ A_k – скорости подачи (1/с) и концентрации (мольн. доли) реагентов на входе и выходе реактора. Если в реакции выполняются линейные стехиометрические законы сохранения (ЛСЗС) вида

где α_jk и C_j – константы, то размерность системы (2) можно уменьшить на число этих ЛСЗС. Как показано в работе [40], в соответствии с правилом стехиометрии Гиббса [41] точное число ЛСЗС (зависимых реагентов) равно

где R_s ≡ rank(a_+ik – a_–ik) – ранг матрицы стехиометрических коэффициентов. Выразим c помощью ЛСЗС N_s зависимых реагентов через остальные (независимые) и исключим их из выражений (2), (3). С учетом этого далее будем считать, что уравнения (2) включают только независимые реагенты, число которых равно K. Отметим, что в качестве независимых могут быть выбраны любые, наиболее удобные для измерений реагенты.

В закрытых системах стационарные состояния (с.с.) определяются только константами скоростей стадий. В открытых системах координаты с.c. зависят также от скорости подачи реагентов и начальных условий (н.у.):

Соотношения (5) означают, что в открытых системах, в отличие от закрытых, эксперименты с различными н.у. (мультиэксперименты) характеризуются разными с.с. и могут использоваться в качестве статистической базы для решения обратных задач. Проведем n = 1, 2, …, N мультиэкспериментов и измерим координаты с.с. в каждом из них

Подставим (6) в (2) и получим систему K × N линейных по k_±i уравнений:

где k = 1, …, K, n = 1, 2, …, N. Эти уравнения разрешимы только тогда, когда

При K × N = 2s и Δ ≠ 0 система (7) не вырождена и константы скоростей всех стадий определяются однозначно:

где Δ = Δ(a_±ik, $q_{n}^{0},$ $A_{{kn}}^{0},$ A_kn) и Δ_±i = Δ_±i(a_±ik, $q_{n}^{0},$ $A_{{kn}}^{0},$ A_kn) – главный и вспомогательный определители системы (7). При этом некоторые из них могут оказаться не физичными, а условия физичности имеют вид

При K × N < 2s система (7) вырождена и константы скоростей стадий определяются неоднозначно:

где Δ_KN = Δ(a_±ik, $q_{n}^{0},$ $A_{{kn}}^{0},$ A_kn, k_±i*) ≠ 0 и Δ_±i, KN = = Δ_±i(a_±ik, $q_{n}^{0},$ $A_{{kn}}^{0},$ A_kn, k_±i*) – главный и вспомогательный определители системы из K × N уравнений; k_±i*, i* = 1, …, 2s – K × N – независимые константы (их значения не определяются однозначно и могут быть заданы произвольно, например, k_±i* = 0), условия физичности которых аналогичны (10). При K × N > 2s или невыполнимости хотя бы одного из условий (10) обратная задача не имеет физичных решений при данном выборе н.у. для мультиэкспериментов. В этих случаях необходимо выбрать другой набор экспериментальных точек.

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

РЕЗУЛЬТАТЫ И ИХ ОБСУЖДЕНИЕ

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

Для этой схемы уравнения (2) для реагентов A, B, C и D соответственно запишутся как

Стехиометрическая матрица схемы (1.1) содержит две строки (по числу стадий) и четыре столбца (по числу реагентов) (1 –1 0 0; 0 1 –1 –1). Ее определитель det(a_ik – a_–ik) ≠ 0, а ранг R_s = 2, и, согласно (4), точное число ЛСЗС (зависимых реагентов) N_s = 4 – 2 = 2. Следовательно, в системе (1.2), (1.3) есть два ЛСЗС, которые нетрудно записать: A + B + C = 1 и A + B + D = 1. Выберем в качестве независимых реагентов, например, А и С, выразим с помощью ЛСЗС остальные и получим систему уравнений для независимых реагентов:

Учитывая, что реакция (1.1) модельная и для нее нет экспериментальных данных, зададим произвольно $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 1, $k_{{ + 2}}^{*}$ = 1, $k_{{ - 2}}^{*}$ = 1 и будем считать их “истинными” значениями констант скоростей стадий реакции. Выберем с учетом ЛСЗС “естественные” начальные условия для первого эксперимента: q = q⁰ = 1, $A_{1}^{0}$ = 1, $B_{1}^{0}$ = 0, $C_{1}^{0}$ = 0, $D_{1}^{0}$ = 0. Выберем также с учетом ЛСЗС н.у. для второго эксперимента: q = q⁰ = 1, $A_{2}^{0}$ = 3/4, $B_{2}^{0}$ = 1/4, $C_{2}^{0}$ = 0, $D_{2}^{0}$ = 0. Решим численно систему (1.4), (1.5) для выбранных н.у., найдем координаты двух с.с. (A₁, С₁) = (0.6065, 0.1805), (A₂, С₂) = = (0.2222, 0.3333) и будем считать их экспериментальными значениями. Подставим все эти значения в (1.4), (1.5) и получим систему линейных уравнений по константам скоростей стадий:

Оценим погрешности метода, решая эту систему с разными н.у. (табл. 1).

Из табл. 1 видно, что найденные значения констант практически не зависят от выбора н.у. Устойчивость метода оценивалась вариацией “истинных” значений констант. Например, для “истинных” значений $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 2, $k_{{ - 2}}^{*}$ = 0.2 метод дает k₊₁ = 1, k_–1 = 0.1, k₊₂ = 2, k_–2 = 0.2 (R = 0.0011). Для более “жестких” значений $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 10, $k_{{ - 2}}^{*}$ = 100 получим k₊₁ = 1, k_–1 = 0.1, k₊₂ = 9.9365, k_–2 = 99.3418 (R = 0.22). Как видно, метод устойчив.

Влияние шума на ошибки вычислений показано в табл. 2. Из этой таблицы видно, что с ростом уровня шума до 20% погрешность определения констант не превышает 10%. Соответственно, решениями обратной задачи (с учетом шума) можно считать интервалы k₊₁ ∈ [0.6806, 1], k_–1 ∈ [1, 1 ], k₊₂ ∈ [1.0094, 1.1899], k_–2 ∈ [0.9997, 0.9997], которые близки к “истинным” значениям констант скоростей стадий.

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

Уравнения (2) для этой реакции запишутся в виде

Стехиометрическая матрица схемы (2.1) содержит три строки и три столбца (1 –1 0; 1 0 –2; –1 1 1). Ее ранг R_s = 3 и, согласно выражению (4), точное число ЛСЗС (зависимых реагентов) N_s = 3 – 3 = 0. Следовательно, в этой реакции нет ЛСЗС и все реагенты независимы. Зададим “истинные” значения констант скоростей стадий $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 1, $k_{{ + 2}}^{*}$ = 1, $k_{{ - 2}}^{*}$ = 1, $k_{{ + 3}}^{*}$ = 1, $k_{{ - 3}}^{*}$ = 1 и произвольные н.у. для двух экспериментов, например: q = q⁰ = 1, $A_{1}^{0}$ = 1, $B_{1}^{0}$ = 0, $C_{1}^{0}$ = 0 и $A_{2}^{0}$ = 1, $B_{2}^{0}$ = 0, $C_{2}^{0}$ = 0. Рассчитаем для них координаты двух с.с. (A₁, B₁, С₁) = = (0.4769, 0.3709, 0.5695) и (A₂, B₂, С₂) = (0.2308, 0.6333, 0.3078); примем их за экспериментальные данные и подставим в уравнения (2.2)–(2.4):

Оценим погрешности метода, решая эту систему с разными н.у. (табл. 3).

Из табл. 3 видно, что расчетные значения констант слабо зависят от начальных условий. Для более жестких “истинных” значений констант, например для $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 2, $k_{{ - 2}}^{*}$ = 0.2, $k_{{ + 3}}^{*}$ = 3, $k_{{ - 3}}^{*}$ = 0.3 ($A_{1}^{0}$ = 1, $A_{2}^{0}$ = 3/4), метод дает k₊₁ = = 0.9994, k_–1 = 0.0990, k₊₂ = 2.0007, k_–2 = 0.2004, k₊₃ = 3.0059, k_–3 = 0.3017 (R = 0.0356). Для еще более “жестких” значений $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 0.1, $k_{{ + 2}}^{*}$ = 10, $k_{{ - 2}}^{*}$ = 100, $k_{{ + 3}}^{*}$ = 1000, $k_{{ - 3}}^{*}$ = 0.001 метод теряет устойчивость (решение не физично): k₊₁ = 1, k_–1 = = 0.1056, k₊₂ = 9.9698, k_–2 = 99.6740, k₊₃ = 757.0308, k_–3 = –0.2408 (R = 8.1858). В этом случае следует выбрать другие начальные условия. Например, для н.у. $A_{1}^{0}$ = 1, $A_{2}^{0}$ = 0 и этих же истинных значений констант получим физичное решение k₊₁ = 1, k_–1 = 0.1016, k₊₂ = 10.0035, k_–2 = 100.0366, k₊₃ = = 999.5098, k_–3 = 0.0009 (R = 0.9202).

Влияние шума на ошибки вычислений показано в табл. 4. Из этой таблицы видно, что с ростом уровня шума до 4% погрешность определения констант не превышает 10%. Соответственно, решениями обратной задачи (с учетом шума) можно считать интервалы k₊₁ ∈ [0.85, 1], k_–1 ∈ [1, 1 ], k₊₂ ∈ [0.75, 1], k_–2 ∈ [0.70, 1], k₊₃ ∈ [0.60, 1], k_–3 ∈ ∈ [0.97, 1], которые близки к “истинным” значениям констант скоростей стадий.

Пример 3. Рассмотрим еще более сложную схему реакции:

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

Для схемы (3.1) R_s = 4 и N_s = 0. Следовательно, в этой системе нет зависимых реагентов и все реагенты независимы. Зададим “истинные” значения констант скоростей стадий: $k_{{ + 1}}^{*}$ = 1, $k_{{ - 1}}^{*}$ = 1, $k_{{ + 2}}^{*}$ = 1, $k_{{ - 2}}^{*}$ = 1, $k_{{ + 3}}^{*}$ = 1, $k_{{ - 3}}^{*}$ = 1, $k_{{ + 4}}^{*}$ = 1, $k_{{ - 4}}^{*}$ = 1 и начальные условия для двух экспериментов, например, q = q⁰ = 1, $A_{1}^{0}$ = 1, $B_{1}^{0}$ = 0, $C_{1}^{0}$ = 0, $D_{1}^{0}$ = 0 и $A_{1}^{0}$ = 0, $B_{1}^{0}$ = 1, $C_{1}^{0}$ = 0, $D_{1}^{0}$ = 0. Рассчитаем для них координаты с.с., подставим эти значения в (3.2)–(3.5) и получим восемь линейных уравнений для определения восьми неизвестных констант скоростей стадий. Результаты решения этой системы уравнений приведены в табл. 5. Из этой таблицы видно, что при уровне шума до 0.5% ошибка метода не превышает 8%. При дальнейшем росте уровня шума метод теряет устойчивость. Соответственно, решениями обратной задачи (с учетом шума) можно считать интервалы k₊₁ ∈ [1, 1.0023], k_–1 ∈ ∈ [0.9556, 1.0001], k₊₂ ∈ [0.9999, 1.1239], k_–2 ∈ [0.9998, 1.2661], k₊₃ ∈ [0.8646, 1.0001], k_–3 ∈ [0.9996, 1.0544], k₊₄ ∈ [0.7593, 1.0002], k_–4 ∈ [0.1289, 1.0010], которые близки к “истинным” значениям констант скоростей стадий.

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

где A = [Cp₂ZrH₂ ⋅ ${\text{ClAlBu}}_{{\text{2}}}^{{(i)}}$]₂, B = [Cp₂ZrH₂ ⋅ ⋅ ${\text{ClAlBu}}_{{\text{2}}}^{{(i)}}$], C = ${\text{HAlBu}}_{2}^{{(i)}},$ D = [Cp₂ZrH₂ ⋅ ${\text{HAlBu}}_{2}^{{(i)}} \cdot $ ⋅ ${\text{ClAlBu}}_{{\text{2}}}^{{(i)}}$]₂. Для этой схемы в работе [27] оптимизационными методами найдены “истинные” значения констант скоростей стадий (пересчитаны в 1/с): k₊₁ ≈ 0.066, k_–1 ≈ 0.0829, k₊₂ ≈ 0.263, k_–2 ≈ 0.0287, описывающие кинетику данной реакции со стандартной ошибкой 5%.

Рассчитаем эти константы без применения методов оптимизации с помощью соотношений (2)–(10). Запишем для схемы (4.1) уравнения (2):

Согласно (3), (4), в системе (4.2)–(4.5) есть два независимых ЛСЗС: 2A + B + D = 2, C + D = 1, т.е. два независимых реагента (K = 2). Выберем независимыми реагентами А и D, исключим из (4.2)–(4.5) с помощью ЛСЗС зависимые реагенты и получим эквивалентную систему уравнений с четырьмя неизвестными константами скоростей стадий:

Согласно выражениям (8), (9), для однозначного определения всех неизвестных этой системы необходимое число экспериментов N = s/K = 2. Выберем по данным работы [26] результаты двух экспериментов: A₁ ≈ 0.33, D₁ ≈ 0.82, q₁ = $q_{1}^{0}$ = 0 и A₂ ≈ 0.38, D₂ ≈ 0.76, q₂ = $q_{2}^{0}$ = 0.01 (остальные н.у. одинаковы: A⁰ = 1, D⁰ = 0). Подставим эти значения в (4.6), (4.7) и получим четыре линейных уравнения вида (7) с четырьмя неизвестными константами скоростей стадий:

Решение этой системы, согласно (9), примет вид

Отметим, что для этих экспериментов условия физичности (10) и н.у. выполняются:

Численные значения констант (4.11)–(4.14) при разных уровнях шума приведены в табл. 6. Из этой таблицы видно, что с ростом уровня шума до 20% погрешности метода не превышают 10%. При дальнейшем повышении уровня шума некоторые константы становятся нефизичными. Соответственно, результатом решения обратной задачи являются интервалы k₊₁ ∈ [0.03, 0.06], k_–1 ∈ [0.07, 0.08], k₊₂ ∈ [0.26, 0.51], k_–2 ∈ [0.03, 0.04], которые близки к интервалам их изменений k₊₁ ∈ [0.02, 0.08], k_–1 ∈ ∈ [0.08, 0.35], k₊₂ ∈ [0.07, 0.22], k_–2 ∈ [0.02, 0.08], рассчитанным в работах [26, 27] (пересчитаны в 1/с).

ЗАКЛЮЧЕНИЕ

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

Автор выражает благодарность В.Х. Федотову за полезные обсуждения работы.