A General Approach to Molecular Reconstruction Accuracy Estimation

Full Text

Сокращения и обозначения: x_i – мольная доля i-го компонента; H – информационная энтропия Шеннона; M – молекулярная масса; ω – массовая доля углерода; J – Якобиан (матрица производных расчетных свойств по параметрам модели); A, B, C, D – элементы матрицы Якобиана 2 × 2 в примерах; E – математическое ожидание; N – количество соединений в смеси; Var – дисперсия; p – свойство индивидуального компонента; P – вероятность; σ – среднеквадратичное отклонение; $\overline{X}$ – среднее значение величины X; X^T – транспонированная матрица Х; X’ – производная Х по своему аргументу; X’_p – производная Х по p для функций многих переменных.

Введение

Моделирование кинетики химических реакций давно стало неотъемлемой частью как теоретического познания мира, так и практических приложений.

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

В связи с этим для моделирования, например, реакций в тяжелых нефтяных фракциях вводятся дополнительные предположения, позволяющие сократить число определяемых параметров. Достаточно популярный подход заключается в “агрегировании” разных соединений в групповые компоненты. Тем самым ставится знак равенства между их кинетическими параметрами, что можно рассматривать как достаточно простой способ моделировать эти параметры. Например, в работе [1] именно таким образом решают проблему большого количества компонентов при моделировании каталитического риформинга бензина, а в работе [2] биомассу представляют в виде трехкомпонентной смеси.

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

Идея моделировать значения, традиционно находимые экспериментально, распространяется не только на кинетические параметры, но и на концентрации компонентов, что особенно полезно при применении детальных кинетических моделей [4]. Такие методы в литературе часто называют методами “молекулярной реконструкции” [5], и уже сейчас существует достаточно много разных вариантов, отличающихся предположениями при построении состава, способами определения параметров и тем, что понимается под “компонентом” смеси (начиная от молекул, псевдокомпонентов и узких фракций и заканчивая другими сложными смесями).

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

Теоретическая часть

Обобщенный взгляд на молекулярную реконструкцию

Как отмечалось во введении, ключевым предположением молекулярной реконструкции является существование функциональной связи между (неизвестными) концентрациями и (известными) характеристиками компонентов:

$x_i=\tilde{f}\left(p_i\right)$, (1)

где x_i – концентрация i-го компонента, p_i – свойство i-го компонента, а $\tilde{f}$ – функция.

Что можно сказать о свойствах такой функции? Поскольку сумма всех концентраций должна равняться единице, то

$\tilde{f}\left(p_i\right)=x_i=\frac{x_i}{1}=\frac{x_i}{\sum x_j}=\frac{f\left(p_i\right)}{\sum f\left(p_j\right)}$, (2)

где f – функция, описывающая концентрации с точностью до постоянного множителя (в отличии от $\tilde{f}$). Такая функция должна быть одного знака и не равняться нулю сразу во всех рассматриваемых точках.

Разные методы реконструкции получаются при разных стратегиях выбора f(x). Например, если все концентрации ненулевые, то достаточно естественно искать решение в форме e^g(x), что часто делается в методе реконструкции максимизации энтропии [6]. Стохастическая реконструкция [7, 8] получается, если задать f(x) = F(x) – F(x – 1) для кумулятивного распределения F непрерывной случайной величины (часто гамма-распределения или экспоненциального).

В случае, если в распоряжении имеются данные о нескольких характеристиках, то f(p₁, p₂, …, p_n) часто представляют в виде произведения функций одной переменной, то есть:

$f\left(p_1,p_2,\dots ,p_n\right)=f_1\left(p_1\right)\times f_2\left(p_2\right)\times \dots \times f_n\left(p_n\right)$. (3)

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

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

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

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

В связи с вышесказанным представляется, что для использования реконструированного состава для моделирования кинетики необходимо обеспечить зависимость параметров модели от расчетных свойств.

Предлагаемый метод оценки точности параметров модели

Оценки расчетных свойств находятся однозначно исходя из набора концентраций, а концентрации в свою очередь обычно однозначно определяются на основе параметров модели. Тем самым можно установить в соответствие расчетные свойства смеси и набор параметров модели:

$G\left(u\right)=m$, (4)

где G(u) – функция, ставящая в соответствие параметры модели и расчетные свойства, а m – экспериментально измеренные свойства образца.

Реконструировать образец – значит найти такое u, чтобы G(u) равнялось заданному экспериментальному набору свойств. Стоит отметить, что измерения (m) в реальности заданы с какой-то погрешностью Δm, а значит, нужно ожидать, что и параметры модели (u) могут содержать неопределенность (Δu):

$G\left(u+ \Delta u\right)=m+ \Delta m$. (5)

Если считать, что отклонения невелики, то возможно разложить зависимость расчетных свойств от параметров модели в ряд Тейлора в окрестности u₀:

$G\left(u_0\right)+J\Delta p=m+\Delta m$, (6)

где J – матрица Якоби.

При реконструировании образца обычно находится такие u₀, что G(u₀) = m, а значит JΔu = Δm. Стоит отметить, что Δm – достаточный понятный случайный вектор. В отсутствие систематической ошибки можно ожидать, что каждый элемент этого вектора описывает случайную ошибку с нулевым математическим ожиданием. Более того, в силу независимости измерений различных свойств можно считать, что недиагональные элементы ковариационной матрицы равны нулю. Диагональные элементы являются квадратами среднеквадратичных отклонений и обычно известны из методик.

Используя определение ковариационной матрицы, можно получить связь между ковариационной матрицей параметров модели и ковариационной матрицей измерений:

$\text{cov}\left(\Delta m\right)=E\left(\left(\Delta m-E\left(\Delta m\right)\right){\left(\Delta m-E\left(\Delta m\right)\right)}^T\right)=E\left(\Delta m \Delta m^T\right)$, (7)

$E\left(\Delta m \Delta m^T\right)=E\left(J\Delta u \Delta u^T{J}^T\right)=J\text{cov}\left(\Delta u\right)J^T$. (8)

Для обратимого якобиана возможно получить достаточно простое выражение для ковариационной матрицы параметров:

$\text{cov}\left(\Delta u\right)=J^{-1}\text{cov}\left(\Delta m\right){J^{-1}}^T$. (9)

Если якобиан необратим, то, вероятно, наилучшей оценкой ковариационной матрицы параметров может служить аналогичное выражение с псевдообратной матрицей:

$\text{cov}\left(\Delta u\right)=J^{+}\text{cov}\left(\Delta m\right){J^{+}}^T$. (10)

На диагонали этой матрицы можно найти дисперсии для всех параметров, а недиагональные элементы характеризуют взаимосвязь между параметрами.

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

На практике для вычисления якобиана хорошо подходит прямое автоматическое дифференцирование.

Реконструкция максимизацией энтропии

В терминах вышеописанного подхода реконструкция максимизации энтропии ищет f(p) в виде e^g(p). Функция g(p) описывает все налагаемые ограничения на состав, будь то строгие/нестрогие ограничения на взвешенные средние или дисперсии. При наличие большего количества ограничений f(p), как отмечалось ранее, перемножаются, а значит g(p) просто складываются.

Для строгих ограничений на среднее g(p) = –λp, где p представляет собой вес, с которым данный компонент входит в взвешенное среднее, а λ – это настраиваемый параметр модели (одинаковый для всех компонентов).

Для нестрогих ограничений на среднее g(p) = ε(p/σ), где ε – это настраиваемый параметр модели, и σ – это заданное значение среднеквадратичного отклонения, характеризующего то, насколько мы не уверены в точности экспериментально полученного значения для среднего.

Для строгих ограничений на дисперсию g(p) = –ν(p – p₀)², где ν – параметр модели, а р₀ – это среднее.

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

$H= -\underset{i=1}{\overset{N}{}}{x}_i\ln x_i$. (11)

Задача максимизации решается с учетом налагаемых ограничений, с помощью которых включаются экспериментальные данные. Задачи условной оптимизации могут быть превращены в задачи безусловной оптимизации в пространстве большей размерности с использованием метода множителей Лагранжа. Часто в литературе рассматривают строгие и нестрогие ограничения на среднее и среднеквадратичное отклонения. В терминах оптимизационных задач они записываются следующим образом:

$H= -\underset{i=1}{\overset{N}{}}{x}_i\ln x_i+ \mu \left(1- \underset{i=1}{\overset{N}{}}{x}_i\right)+ \underset{s=1}{\overset{S}{}}{\lambda }_s\left(f_s- \underset{i=1}{\overset{N}{}}{x}_i{f}_{i,s}\right)$, (12)

$H= -\underset{i=1}{\overset{N}{}}{x}_i\ln x_i+ \mu \left(1- \underset{i=1}{\overset{N}{}}{x}_i\right)-\frac{1}{2}\underset{k=1}{\overset{K}{}}\frac{{\left(f_j- {{\displaystyle \sum }}_{i=1}^N{x}_i{f}_{i,j}\right)}^2}{{\sigma }^2}$, (13)

$H= -\underset{i=1}{\overset{N}{}}{x}_i\ln x_i+ \mu \left(1- \underset{i=1}{\overset{N}{}}{x}_i\right)+ \underset{t=1}{\overset{T}{}}{\nu }_r\left({\sigma }_r^2- \underset{i=1}{\overset{N}{}}{x}_i{\left(f_t-f_{i,t}\right)}^2\right)$, (14)

где f с одним индексом – это значение свойства для смеси, f с двумя индексами – это свойства чистых компонентов, S и K – это количества строгих и нестрогих ограничений на взвешенное среднее, T – количество строгих ограничений на среднеквадратичное отклонение, а μ, λ, ν – множители Лагранжа.

Первое слагаемое представляет собой информационную энтропию, второе слагаемое накладывает ограничение на сумму концентраций. Третье слагаемое описывает строгие/нестрогие ограничения на среднее/строгое ограничение на среднеквадратичное отклонение соответственно.

Можно заметить, что в такой форме они все зависят от бо́льшего числа переменных, чем количество искомых концентраций. Однако для задач с такими ограничениями (и любыми комбинациями таковых) существует достаточно простая зависимость концентраций компонентов от множителей Лагранжа x = e^{g(λ, ε, ν)}/Z (Z – нормировочный множитель), не включающая другие концентрации. Поэтому возможно полностью избавиться от концентраций в оптимизационной задаче, заметно сокращая ее размерность. С учетом этого можно записать:

$H=\ln Z+ \underset{s=1}{\overset{S}{}}{\lambda }_s{f}_j$, (15)

$H=\ln Z+\underset{k=1}{\overset{K}{}}\left(\frac{1}{2}{\epsilon }_k^2-\frac{f_k}{{\sigma }_k}{\epsilon }_k\right)$, (16)

$H=\ln Z+ \underset{t=1}{\overset{T}{}}{\nu }_t{\sigma }_t^2$. (17)

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

В общем случае такой “фокус” не работает, и g(p) будет включать в себя концентрации компонентов в явном виде, что может превратить задачу в достаточно нетривиальную.

Рассмотрим общий случай, при котором оценка свойства определяется функцией L(x):

$H= -\underset{i=1}{\overset{N}{}}{x}_i\ln x_i+ \mu \left(1- \underset{i=1}{\overset{N}{}}{x}_i\right)+ \underset{s=1}{\overset{S}{}}{\lambda }_s\left(f_s- L\left(x\right)\right)$. (18)

Частная производная этого выражения по концентрации i-го компонента:

$\frac{\partial H}{\partial x_i}= -\ln x_i-1-\mu - \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}$. (19)

В условиях оптимума производная будет равна нулю, и тогда:

$-\ln x_i-1-\mu - \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}=0$, (20)

$\ln x_i=-\mu -1- \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}$, (21)

$x_i={\text{e}}^{-\mu -1- \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}}$, (22)

$x_i{\text{e}}^{\mu +1}={\text{e}}^{- \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}}$. (23)

Просуммировав такие выражения для каждого i:

${\text{e}}^{\mu +1}\underset{i=1}{\overset{N}{}}{x}_i=\underset{i=1}{\overset{N}{}}{\text{e}}^{- \underset{s=1}{\overset{S}{}}{\lambda }_s\frac{\partial L}{\partial x_i}}$. (24)

Поскольку сумма концентраций равна 1, то уравнение (22) можно переписать в виде:

$x_i=\frac{{\text{e}}^{- {{\displaystyle \sum }}_{s=1}^S{\lambda }_s\frac{\partial L}{\partial x_i}}}{{\text{e}}^{\mu +1}}= \frac{{\text{e}}^{- {{\displaystyle \sum }}_{s=1}^S{\lambda }_s\frac{\partial L}{\partial x_i}}}{{\displaystyle {\sum }_{i=1}^N{\text{e}}^{- {\sum }_{s=1}^S{\lambda }_s\frac{\partial L}{\partial x_i}}}}$. (25)

В наших прошлых обозначениях это просто значит, что в общем случае для строго ограниченного среднего g(x) = –λL’(x), где L’ – это производная правил смешения по концентрации компонента. Если правила смешения сложнее, чем линейные, то производная по концентрации не будет константой, что в конечном счете и позволяло упрощать задачу в рассмотренных выше случаях.

Стохастическая реконструкция

Другой достаточно популярный подход к моделированию составов представляет компоненты как комбинацию “структурных атрибутов”, которые, как понятно из названия, характеризуют структуру молекул. Концентрации компонентов определяются как функции этих самых структурных атрибутов. Структурные атрибуты могут принимать как числовые значения (длина цепи парафина, количество/длина заместителей, количество колец и т.д.) так и категориальные (тип молекулы).

В качестве основы зависимости концентрации от структурных атрибутов используется кумулятивные распределения случайных величин. Например, можно задавать f(x) как

$f\left(x\right)=F\left(x\right)-F\left(x-1\right)$ (26)

для произвольного кумулятивного распределения F.

Например, если структурный атрибут может принимать значения от 1 до Q, то вероятность получить q будет:

$P\left(X=q\right)= \frac{F\left(q\right) - F\left(q-1\right)}{F\left(Q\right)-F\left(0\right)}$. (26)

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

В классической реализации стохастической реконструкции составы получаются именно по такой схеме с генерацией случайных чисел. Концентрации определяются как наблюдаемые частоты (т.е. после генерации N “молекул” таким способом, каждой присваивается концентрация 1/N).

Экспериментальная часть

Метод максимизации энтропии

Предположим, что имеется смесь нормальных парафинов (для ясности, скажем, С₁₀–С₂₀) и известна средняя молекулярная масса этой смеси. Опираясь на эти данные, почти наверняка невозможно определить состав точно, но можно его оценить с помощью метода максимизации энтропии.

Средняя молекулярная масса смеси представляет собой взвешенную сумму молекулярных масс компонентов. Используя строгие ограничения f(M) = e^–λM и варьируя значение λ, можно добиться почти любого значения средней молекулярной массы.

Используя рис. 1, можно понять, что погрешность в молекулярной массе в разной степени влияет на погрешность в описании параметров. Достаточно широкий диапазон значений молекулярных масс может отвечать очень узкому диапазону параметров и, наоборот, погрешность в определении свойств может соответствовать очень широкому диапазону параметров модели и, как следствие, концентраций компонентов. Так, в вышеприведенном примере, если мы знаем, что средняя молекулярная масса находится между 270 и 280, оценка концентрации С₂₀ в смеси будет варьироваться от 53 до 88%.

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

Несложно представить ситуацию, в которой скорости реакций тяжелых парафинов будут значительно выше, чем скорости легких. Соответственно, на результаты моделирования сильно повлияет то обстоятельство, составляют ли тяжелые парафины 90% всей смеси или 50%.

Предлагаемый метод оценки неопределенности параметров в данном случае утверждает:

${\sigma }_u^2= \frac{1}{M^{\text{'}}\left(\lambda \right)}{\sigma }_M^2\frac{1}{M^{\text{'}}\left(\lambda \right)}= {\left(\frac{{\sigma }_M}{M^{\text{'}}\left(\lambda \right)}\right)}^2$, (27)

${\sigma }_u=\frac{{\sigma }_M}{\left|M\text{'}\left(\lambda \right)\right|}$. (28)

Для удобства вычислений можно показать, что производная при таких правилах смешения равна дисперсии с точностью до знака:

$x_i=\frac{{\text{e}}^{-\lambda M_i}}{\sum {\text{e}}^{-\lambda M}}$, (29)

${x^{\text{'}}}_i=\frac{-M_i{\text{e}}^{-\lambda M_i} \sum {\text{e}}^{-\lambda M}+ {\text{e}}^{-\lambda M_i}\sum M{\text{e}}^{-\lambda M}}{{\left(\sum {\text{e}}^{-\lambda M}\right)}^2}= -M_i{x}_i+x_i\overline{M}$, (30)

${\left(\underset{i=1}{\overset{N}{}}{M}_i{x}_i\right)}^{\text{'}}= -\underset{i=1}{\overset{N}{}}{M}_i^2{x}_i+ {\overline{M}}^2= {\overline{M}}^2-\overline{M^2}= -\text{Var}\left(M\right)$. (31)

Можно построить зависимость среднеквадратичного отклонения параметра (σ) от значения параметра (рис. 2), а также среднюю абсолютную ошибку в определении концентраций (рис. 3).

Рис. 2. Зависимость погрешности находимого параметра от значения параметра.

Рис. 3. Зависимость средней абсолютной ошибки определения концентрации от параметра модели.

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

При использовании нестрогих ограничений (f(x) = e^εx/σ) состав становится ближе к эквимолярному по сравнению с теми же значениями для строгих ограничений (рис. 4). Можно заметить, что показатель экспоненты в этой модели не содержит –1 из строгих ограничений. Это не имеет принципиального эффекта на решение.

Рис. 4. Зависимость расчетной молекулярной массы смеси в зависимости от значения параметра модели для точных и неточных ограничений.

Несложно заметить, что четкие ограничения с параметром λ соответствуют нечетким ограничениями с –ε/σ (т.е. отличаются на постоянный множитель), а значит, производная молекулярной массы по ε будет:

$M^{\text{'}}\left(\epsilon \right)= \frac{\text{d}M}{\text{d}\lambda }\frac{\text{d}\lambda }{\text{d}\epsilon }=\frac{\text{Var}\left(M\right)}{\sigma }$. (32)

Тем самым,

${\sigma }_u^2=\frac{1}{M^{\text{'}}\left(\epsilon \right)}{\sigma }_M^2\frac{1}{M^{\text{'}}\left(\epsilon \right)}=\frac{{\sigma }_M^2 {\sigma }^2}{\text{Var}{\left(M\right)}^2}$, (32)

${\sigma }_u=\frac{\sigma {\sigma }_M}{\text{Var}\left(M\right)}$. (33)

В сущности это означает, что неопределенности перемножаются.

Введем в рассмотрение дополнительные измерения, скажем, среднюю долю углерода в смеси. Обычно, f(M, ω) = f₁(M) × f₂(ω) = e^{-λ1M – λ2ω}. Вычислим Якобиан. Для начала определимся с производными концентраций по параметрам модели:

$\frac{\partial x_i}{\partial {\lambda }_1}= {{\left(\frac{{\text{e}}^{-{\lambda }_1{M}_i- {\lambda }_2{\omega }_i}}{\sum {\text{e}}^{-{\lambda }_1{M}_i- {\lambda }_2{\omega }_i}}\right)}^{\text{'}}}_{{\lambda }_1}= -M_i{x}_i+x_i\overline{M}=x_i\left(\overline{M}- M_i\right)$, (34)

$\frac{\partial x_i}{\partial {\lambda }_2}= {{\left(\frac{{\text{e}}^{-{\lambda }_1{M}_i- {\lambda }_2{\omega }_i}}{\sum {\text{e}}^{-{\lambda }_1{M}_i- {\lambda }_2{\omega }_i}}\right)}^{\text{'}}}_{{\lambda }_2}= -{\omega }_i{x}_i+x_i\overline{\omega }=x_i\left(\overline{\omega }- {\omega }_i\right)$. (35)

Тогда элементы Якобиана можно найти как

$\frac{\partial \overline{M} }{\partial {\lambda }_1}= {\left(\underset{i=1}{\overset{N}{}}{M}_i{x}_i\right)}^{\text{'}}=-\text{Var}\left(M\right)$, (36)

$\frac{\partial \overline{M} }{\partial {\lambda }_2}={\left(\underset{i=1}{\overset{N}{}}{M}_i{x}_i\right)}^{\text{'}}=\underset{i=1}{\overset{N}{}}{M}_i{x}_i=\underset{i=1}{\overset{N}{}}{M}_i{x}_i\left(\overline{\omega }- {\omega }_i\right)= \overline{M}\left(\overline{\omega }-\omega \right)$, (37)

$\frac{\partial \omega }{\partial \lambda }={\left(\frac{\sum M_i{\omega }_i{x}_i}{\sum M_i{x}_i}\right)}^{\text{'}}= {\left(\frac{\overline{M\omega }}{\overline{M}}\right)}^{\text{'}}=\frac{{\overline{M\omega }}^{\text{'}}\overline{M}-\overline{M\omega }{\overline{M}}^{\text{'}} }{{\overline{M}}^2}=\frac{{\overline{M\omega }}^{\text{'}}- \omega {\overline{M}}^{\text{'}}}{\overline{M}}$, (38)

$\frac{\partial \omega }{\partial {\lambda }_1}=\frac{{\overline{M\omega }}^{\text{'}}+ \omega \text{Var}\left(M\right)}{\overline{M}}=\frac{\omega \overline{M^2}-\overline{M^2\omega } }{\overline{M}}$, (39)

$\frac{\partial \omega }{\partial {\lambda }_2}=\frac{{\overline{M\omega }}^{\text{'}}+ \omega \text{Var}\left(M\right)}{\overline{M}}= \omega \left(\overline{\omega }+ \frac{\text{Var}\left(M\right)- \overline{M{\omega }^2}}{\overline{M}}\right)$. (40)

Тогда Якобиан будет равен

$J=\left[\begin{array}{cc}-\text{Var}\left(M\right)& \overline{M}\left(\overline{\omega }-\omega \right)\\ \frac{\omega \overline{M^2}-\overline{M^2\omega } }{\overline{M}}& \omega \left(\overline{\omega }+ \frac{\text{Var}\left(M\right)- \overline{M{\omega }^2}}{\overline{M}}\right)\end{array}\right]=\left[\begin{array}{cc}A& B\\ C& D\end{array}\right]$, (41)

а ковариационная матрица параметров

$\text{cov}\left(\Delta p\right)=\frac{1}{{\left(\text{det}J\right)}^2}\left[\begin{array}{cc}D& -B\\ -C& A\end{array}\right]\left[\begin{array}{cc}{\sigma }_M^2& 0\\ 0& {\sigma }_{\omega }^2\end{array}\right]\left[\begin{array}{cc}D& -C\\ -B& A\end{array}\right]=\frac{1}{{\left(\text{det}J\right)}^2}\left[\begin{array}{cc}{D}^2{\sigma }_M^2+B^2{\sigma }_{\omega }^2& -\left(DC{\sigma }_M^2+AB{\sigma }_{\omega }^2\right)\\ -\left(DC{\sigma }_M^2+AB{\sigma }_{\omega }^2\right)& C^2{\sigma }_M^2+A^2{\sigma }_{\omega }^2\end{array}\right]$. (42)

А значит,

${\sigma }_{{\lambda }_1}^2=\frac{D^2{\sigma }_M^2+B^2{\sigma }_{\omega }^2}{{\left(\text{det}J\right)}^2}$, (43)

${\sigma }_{{\lambda }_2}^2=\frac{C^2{\sigma }_M^2+A^2{\sigma }_{\omega }^2}{{\left(\det J\right)}^2}$. (44)

Видно, что дисперсия параметров может быть выражена как взвешенная сумма дисперсий измерений.

Может показаться достаточно контринтуитивно, что в выражении для неопределенности в параметре, отвечающем за молекулярную массу, фигурирует произведение неопределенности в экспериментальном определении молекулярной массы и производной массовой доли углерода по параметру, отвечающему за массовую долю углерода, но если переписать это выражение в следующем виде:

${\sigma }_{{\lambda }_1}^2=\frac{{\sigma }_M^2}{{\left(\frac{\det J}{D}\right)}^2}+\frac{{\sigma }_{\omega }^2}{{\left(\frac{\det J}{B}\right)}^2}= \frac{{\sigma }_M^2}{{\left(A-\frac{BC}{D}\right)}^2}+\frac{{\sigma }_{\omega }^2}{{\left(\frac{AD}{B}-C\right)}^2}$, (45)

то можно заметить, что это достаточно естественное обобщение случая с одной переменной.

Ковариация параметров равна:

$\mathrm{cov}\left({\lambda }_1,{\lambda }_2\right)=\frac{-\left(DC{\sigma }_M^2+AB{\sigma }_{\omega }^2\right)}{{\left(\det J\right)}^2}$, (46)

и коэффициент корреляции между параметрами будет равен:

$r\left({\lambda }_1,{\lambda }_2\right)=\frac{-\left(DC{\sigma }_M^2+AB{\sigma }_{\omega }^2\right)}{\sqrt{D^2{\sigma }_M^2+B^2{\sigma }_{\omega }^2}\sqrt{C^2{\sigma }_M^2+A^2{\sigma }_{\omega }^2}}$. (47)

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

Стохастическая реконструкция

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

Предполагая, что длина парафинов подчиняется гамма-распределению, можно добиться любых допустимых значений молекулярной массы и доли углерода в смеси, варьируя среднее и дисперсию (рис. 5).

Рис. 5. Средняя молекулярная масса (а) и массовая доля углерода (б) смеси в зависимости от параметров модели.

Можно заметить достаточно общую тенденцию для обоих свойств: они довольно сильно меняются с изменением среднего и почти не меняются с изменением дисперсии, что не кажется удивительным, поскольку оба свойства в некотором смысле измеряют “среднее” поведение смеси.

В связи с этим сложно ожидать хоть какой-либо точности в определении параметров. Определитель матрицы Якоби в данном случае для большинства значений составляет порядка 10^–15 (т.е. близок к машинной точности), и коэффициент корреляции равен 1 (самое меньшее значение отстояло от единицы не более чем на 5 × 10^–16). Расчетные значения логарифмов среднеквадратичных отклонений при σ_м = 1 г/моль и σ_ω = 0.1 мас. % показаны на рис. 6. Можно отчетливо видеть численную неустойчивость, вызванную почти сингулярным якобианом.

Рис. 6. Зависимость логарифма среднеквадратичного отклонения среднего (а) и дисперсии (б) в зависимости от истинных параметров модели.

Если же использовать другое свойство, которое в большей степени характеризует дисперсию, то ошибка определения параметров будет заметно ниже. Например, если измерить долю отгона при двух температурах, можно найти произведение долей самой легкой и самой тяжелой частей. Несложно заметить, что для того, чтобы такая величина была большой, необходимо значительное количество легкой и тяжелой частей одновременно. Расчеты точности параметров модели в этом случае приведены на рис. 7. Вторым свойством выступало произведение суммарных концентраций С₁₀–С₁₂ и С₁₈–С₂₀.

Результаты и их обсуждение

Сравнивания рис. 6 и 7, можно заметить, что рис. 7 существенной численной неустойчивости не наблюдается. Кроме того, значения оценки погрешностей меньше почти на 16 порядков.

Рис. 7. Зависимость логарифма среднеквадратичного отклонения среднего (a), дисперсии (б) и коэффициент корреляции (в).

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

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

Заключение

Использование реконструированных составов для описания кинетики может очень сильно помочь в построении детальных кинетических моделей и дальнейшей оптимизации процессов, однако необдуманное применение этих подходов способно привести к абсурдным результатам. В связи с этим, необходимо тщательное исследование привлекаемых моделей.

В работе предложен достаточно общий взгляд на методы реконструкции и предложена формула, позволяющая оценивать точность находимых параметров. Несмотря на изобилие математических выкладок в самой работе, применение выдвинутой формулы на практике не должно вызывать трудностей при применении прямого автоматического дифференцирования /(многомерных) дуальных чисел для нахождения Якобиана. Численная нестабильность, возникающая при высокой корреляции выбранных методов анализа, сама по себе может свидетельствовать о малой надежности находимых параметров.

Знание точности параметров позволяет определить, в каких диапазонах могут варьироваться концентрации компонентов, не противореча экспериментальным данным, и соответственно решать, насколько выбранная модель может быть использована для описания кинетики.

Финансирование

Работа выполнена при поддержке Министерства науки и высшего образования Российской Федерации в рамках государственного задания для Института катализа им. Борескова (проект FWUR-2024-0037).

Конфликт интересов

Авторы заявляют об отсутствии конфликта интересов, требующего раскрытия в данной статье.

About the authors

N. A. Glazov

Author for correspondence.
Email: glazov@catalysis.ru
Russian Federation, 5 Lavrentiev Ave., Novosibirsk, 630090

A. N. Zagoruiko

Email: glazov@catalysis.ru
Russian Federation, 5 Lavrentiev Ave., Novosibirsk, 630090

References

Supplementary files

Supplementary Files

Action

1. JATS XML

Download

2. Fig. 1. Dependence of the calculated average molecular mass on the model parameter: a wide range of molecular masses corresponds to a narrow range of parameters (a) and a narrow range of molecular masses corresponds to a wide range of parameters (b).

Download (20KB)

Indexing metadata

3. Fig. 2. Dependence of the error of the found parameter on the parameter value.

Download (12KB)

Indexing metadata

4. Fig. 3. Dependence of the average absolute error in determining concentration on the model parameter.

Download (12KB)

Indexing metadata

5. Fig. 4. Dependence of the calculated molecular mass of the mixture on the value of the model parameter for exact and imprecise constraints.

Download (19KB)

Indexing metadata

6. Fig. 5. Average molecular weight (a) and mass fraction of carbon (b) of the mixture depending on the model parameters.

Download (29KB)

Indexing metadata

7. Fig. 6. Dependence of the logarithm of the standard deviation of the mean (a) and variance (b) on the true parameters of the model.

Download (36KB)

Indexing metadata

8. Fig. 7. Dependence of the logarithm of the standard deviation of the mean (a), variance (b) and correlation coefficient (c).

Download (28KB)

Indexing metadata