Научная статья на тему 'АНАЛИЗ ИДЕНТИФИЦИРУЕМОСТИ ПАРАМЕТРОВ КИНЕТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКИХ РЕАКЦИЙ'

АНАЛИЗ ИДЕНТИФИЦИРУЕМОСТИ ПАРАМЕТРОВ КИНЕТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКИХ РЕАКЦИЙ Текст научной статьи по специальности «Математика»

CC BY
39
9
i Надоели баннеры? Вы всегда можете отключить рекламу.
Ключевые слова
АНАЛИЗ ЧУВСТВИТЕЛЬНОСТИ / АНАЛИЗ ИДЕНТИФИЦИРУЕМОСТИ / КИНЕТИЧЕСКАЯ МОДЕЛЬ

Аннотация научной статьи по математике, автор научной работы — Сафиуллина Л.Ф., Губайдуллин И.М.

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

i Надоели баннеры? Вы всегда можете отключить рекламу.

Похожие темы научных работ по математике , автор научной работы — Сафиуллина Л.Ф., Губайдуллин И.М.

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.
i Надоели баннеры? Вы всегда можете отключить рекламу.

ANALYSIS OF THE IDENTIFICABILITY OF KINETIC MODEL PARAMETERS OF CHEMICAL REACTIONS

Kinetic models of chemical reactions usually consist of ordinary differential equations with many unknown parameters. Some of these parameters are often practically unidentifiable, that is, their values cannot be unambiguously determined from the available data. Therefore, before defining model parameter values, it is important to characterize the subset of identified parameters and their interactions. The ability to successfully estimate model parameters depends on the structure of the model, the number of parameters, and the set of experimental data. A set of parameters which evaluation is impossible may be adjusted by applying nominal values, or the model may be modified in order to exclude those parameters at all. Alternatively, the developer of the model may propose additional experiments that may increase the number of identified parameters. The authors of the article present methods to determine a subset of identifiable model parameters. Practical application of the method is shown on the example of a model reaction of alpha-pinene isomerization. Using the DAISY software package for the analysis of the structural identifiability of the model, the authors made sure that the mathematical model is identifiable. The sensitivity analysis of all five parameters of the mathematical model was carried out by the orthogonal method and the eigenvalue method. According to the results of sensitivity and identifiability analyses, the values of model parameters satisfactory correspond to the values of experimental data.

Текст научной работы на тему «АНАЛИЗ ИДЕНТИФИЦИРУЕМОСТИ ПАРАМЕТРОВ КИНЕТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКИХ РЕАКЦИЙ»

УДК 519.6+544.4

DOI: 10.33184/bulletin-bsu-2022.2.8

АНАЛИЗ ИДЕНТИФИЦИРУЕМОСТИ ПАРАМЕТРОВ КИНЕТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКИХ РЕАКЦИЙ

© Л. Ф. Сафиуллина1*, И. М. Губайдуллин2

1Башкирский государственный университет Россия, Республика Башкортостан, 450076 г. Уфа, ул. Заки Валиди, 32.

2Институт нефтехимии и катализа РАН Россия, Республика Башкортостан, 450075 г. Уфа, пр. Октября, 141.

*ЕтаИ: safiullinalf@gmail.com

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

Ключевые слова: анализ чувствительности, анализ идентифицируемости, кинетическая модель.

Введение

Анализ идентифицируемости позволяет установить возможность определения значений неизвестных параметров модели [1-3]. Выделяют структурный и практический анализ идентифицируемости. Структурный анализ идентифицируемости решает, являются ли параметры модели однозначно определяемыми на основе структуры модели [4]. Параметр модели к является структурно идентифицируемым, если у(к)=у(к') <=>к=к\ где у - прогнозируемое значение модели. Параметр к является структурно локально идентифицируемым, если для почти любого значения к* существует окрестность V(k*), в которой выполняется указанное выше соотношение. Параметр является глобально идентифицируемым, если отношение выполняется во всем диапазоне значений параметра. Если есть некоторая область с ненулевой мерой, где соотношение не выполняется, параметр к структурно неидентифицируем. Анализ структурной идентифицируемости обычно связан с большими вычислительными затратами, что затрудняет его применение к большим моделям [1]. Кроме того, структурная идентифицируемость является лишь необходимым, но не достаточным условием идентифицируемости. Очень часто структурно идентифицируемый параметр практически неидентифицируем, т.е. его значение не может быть точно определено из-за ограниченности экспериментальных данных. С помощью практического анализа идентифицируемости выделяют параметры, которые не оказывают влияние на экспериментальные наблюдения, а также определяют взаимозависимые параметры. Очевидно, что если параметр не влияет на наблюдаемые вещества, определить его значение невозможно. Вторая ситуация, когда влияние изменения одного параметра на наблюдаемые

вещества может быть компенсировано изменением других параметров, что также может препятствовать идентификации параметра [5].

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

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

Математическая модель

Мы рассматриваем кинетические модели химических реакций, которые могут быть описаны нелинейными обыкновенными дифференциальными уравнениями в следующем виде [6]:

dxi dt

■fr (xt (t, к), к),

(1)

XI(0) = х0,1 = 1,..., п.

Здесь хь - концентрации веществ, fi описывает взаимодействия между переменными состояния (функции скорости реакции). Вектор параметров к содержит (положительные) параметры, например, коэффициенты скорости реакции. Их значения часто неизвестны и должны оцениваться по данным. Переменные модели xi сопоставляются с измеримыми выходными переменными у, также известными как наблюдаемые вещества реакции или предсказания модели. Обозначим через у,у предсказание модели

для /-й наблюдаемой величины в момент времени ]. Соответствующие измеренные данные обозначаются уц.

Определение параметров модели

Целью оценки параметров является определение значений неизвестного вектора параметров к. Обычно это делается путем минимизации расстояния между предсказанием модели у] и измеренными данными уц. Одним из самых часто применяемых вариантов является минимизация средневзвешенной суммы квадратов отклонений [7]:

^ = И^пц (Уц-Уц)2, (2)

где N - количество экспериментов; М - количество наблюдаемых соединений, веса обозначены .

Анализ чувствительности

Анализ чувствительности [8-11] для оценки идентифицируемости неизвестных параметров модели (1), как правило, используется до численного решения задачи идентификации параметров. Исследование чувствительности для математической модели проводится относительно некоторых номинальных параметров, значения которых берутся из литературных источников или доступной статистической информации.

Коэффициенты чувствительности для переменной по параметру к] определяются как предел отношения [12]:

^ (0 = ^ (3)

Уравнения для вычисления получим с помощью частного дифференцирования обеих частей системы уравнения (1) и начальных условий по параметру кр

(4)

= У + 211

М Ьа = 1дх„1) дк{

Здесь ^ - якобиан системы ОДУ, а Ц- - производная правой части по рассматриваемым параметрам.

В качестве начальных условий для этой системы линейных обыкновенных дифференциальных уравнений с постоянными коэффициентами имеем

1 - (0) = 0, {Ф

дХ] (0)

(0) = 1.

(5)

Таким образом, вычисление функций чувствительности для всего комплекса параметров системы (1) сводится к совместному интегрированию уравнений вида (4) для каждого параметра к] (]=1,...,ш) и модели объекта (1).

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

% (0 = '-ЩЧ = % м £ (6)

где Бц - нормированные коэффициенты чувствительности.

Методы анализа идентифицируемости основаны на анализе матрицы, состоящей из нормированных коэффициентов чувствительности [1; 3]:

/811(11) ... Б1т(11-\

5 =

^и(^)

$пт (^1) )

3пт ~)/

(7)

Методы анализа идентифицируемости

1. Ортогональный метод

Ортогональный метод, разработанный Yao [3], хорошо себя зарекомендовал для анализа моделей с большим числом параметров. Данный метод требует заранее заданных значений параметров (номинальные значения либо приближенные значения) и количество и расположение временных точек измерений. Его основная идея состоит в том, чтобы исследовать линейные зависимости столбцов матрицы чувствительности & Таким образом, одновременно можно оценить как чувствительность параметров ко входным данным, так и зависимость между параметрами.

Ниже приведена последовательность шагов метода:

1. Вычислить сумму квадратов элементов для каждого столбца матрицы &

2. Выбрать параметр, сумма квадратов элементов которого имеет наибольшую величину, в качестве первого оцениваемого параметра.

3. Отметить соответствующий столбец как Хь (Ь= 1 для первой итерации).

4. Вычислить §ь = Хь(Х{ХЬ)-1Х[Б.

5. Вычислить перпендикуляры Иь = 5 — §ь.

6. Вычислите сумму квадратов для каждого столбца матрицы Яь. Столбец с наибольшей суммой соответствует следующему идентифицируемому параметру.

7. Выбрать соответствующий столбец в матрице & и увеличить матрицу Хь, включив новый столбец. Обозначим расширенную матрицу Хь +1.

8. Увеличить счетчик итераций на единицу и повторять шаги с 4 по 7 до полного упорядочивания всех параметров модели.

2. Метод собственных значений

Метод собственных значений был впервые предложен Vajda [13], а затем получил дальнейшее развитие [1]. Этот метод основан на свойствах собственных значений и собственных векторов матрицы 5Г5, где S - матрица чувствительности, определенная в (7). Обозначим Х1<Х2< ••• < собственные значения матрицы Н = 5ТБ в неубывающем порядке, а (у1, у2,..., уч) - соответствующие нормированные собственные векторы. Критерий выбора для неидентифицируемых параметров определяется выражением

т = arg(max |ум|).

1<h<q J

(8)

3. Метод корреляционной матрицы

Информацию о корреляции между оцениваемыми параметрами можно получить из ковариационной матрицы. Корреляционная матрица, элементами которой являются приближенные коэффициенты корреляции между i-м и]-м параметром, определяется как:

D — Ч

Mij V^y,

где коэффициенты Cj - это элементы матрицы C:

(9)

С =

(10)

Математическая модель изомеризации альфа-пинена

В качестве тестового примера рассмотрим кинетическую задачу изомеризации альфа-пине-на [14]. Мы хотим оценить 5 констант скорости (к1,...,к5) сложной биохимической реакции, первоначально изученной Боксом и его коллегами [15]. Ниже представлена предполагаемая схема этой гомогенной химической реакции, описывающая термическую изомеризацию а-пинена (х1) в дипентен (Х2) и аллоцимен (хз), которая, в свою очередь, дает а- и р-пиронен (Х4) и димер (Х5).

/с 2 ^3 ^5

А^ В, А^ С,С ^ О, С^ Е, Е^ С, где А - пинен, В - дипентен, С - алло-оцимен , Б -пиролен, Е - димер.

Хантер и др. [14] предположили кинетику первого порядка и вывели следующие линейные уравнения для пяти веществ: йх1

—— = (—к1 — к 2) Х1,

dt

(11)

йх2

И = к1Хъ

йх3

= к 2X1 + к5х5 — к4х3 — к3х3, йх4

И = ^

йх5

= к4х3 — к5х5.

X; (0) = х°,1 = 1, ...,5.

Здесь х1 - А, х2 - В, х3 - С, х4 - Б, х5 - Е.

Бокс и его коллеги [15] обнаружили, что в данных существуют зависимости, поэтому локальные методы оптимизации, а также некоторые глобальные методы попадали в локальную ловушку. После многочисленных исследований данной реакции были определены константы данной реакции [16]: к1=5.936 • 10-5, к2=2.937 • 10-5, кз=1.978 • 10-5, к4=3.084 • 10-4, к5=5.146 • 10-5.

На рис. 1 представлены расчетные кривые, полученные решением системы (11) с известными значениями констант, и экспериментальные данные.

Используя программный пакет DAISY [17] по анализу структурной идентифицируемости модели (11), мы убедились, что математическая модель является идентифицируемой. Далее был проведен анализ чувствительности всех 5-ти параметров математической модели (11) ортогональным методом и методом собственных значений (рис. 2). Набор параметров, чувствительных к ошибкам измерений, упорядоченный в порядке усиления чувствительности, который был получен с помощью обоих методов анализа идентифицируемости, совпадает и является следующим: к2, ки к4, кз, кз. Отметим полное совпадение последовательности чувствительных параметров полученных обоими методами. Для последнего столбца (соответствует параметру кз) значение нормы перпендикуляра составило 0.97. Это достаточно высокое значение, согласно исследованиям [1; 3], значит, все параметры модели являются идентифицируемыми.

10000 20000 зоооо

Время, сек

40000

Рис. 1. Расчетные и экспериментальные данные кинетики изомеризации альфа-пинена.

Для метода корреляционной матрицы положим, что измерения берутся с ошибкой 5 %. Тогда получаем корреляционную матрицу вида (12):

iНе можете найти то, что вам нужно? Попробуйте сервис подбора литературы.

К t

к2 кз

к 4 fc5

0.09 -0.12 0.02 \ 0.03

0.09 1

0.24 0.08 0.15

0.02 0.08 0.26 1

0.82

0.03 0.15 -0.18 0.82 1

(12)

—0.12 0.24 1

0.26 —0.18

Заметим, что коэффициенты корреляции для параметров (к4; к5) равны 0.82 и близки к единице. Это означает, что эти параметры зависят друг от друга и неразличимы для модели (11). Все остальные коэффициенты имеют невысокую корреляцию. Следовательно, для дальнейшего определения параметров данной модели желательно зафиксировать наиболее чувствительные параметры (к4; кз) и определять оставшиеся три. Этот факт также отражает наличие локальных минимумов и проблему определения параметров модели, с которыми сталкивались исследователи.

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

Заключение

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

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 19-37-60003 60003 и в рамках Государственного задания Института нефтехимии и катализа УФИЦ РАН (код (шифр) научной темы - FMRS-2022-0078).

ЛИТЕРАТУРА

1. Miao H., Xia X., Perelson A. S., We H. On Identifiability of nonlinear ode models and applications inviral dynamics // SIAM Rec. Soc. Ind. Appl. Math. 2011. Vol. 53. №.1. P. 3-39.

2. Saccomani M. P., Thomaseth K. The Union between Structural and Practical Identifiability Makes Strength in Reducing Oncological Model Complexity: A Case Study // Complexity. 2018. С. 1-10.

3. Yao K. Z., Shaw B. M., Kou B., McAuley K. B., Bacon D. W. Modeling ethylene/butane copolymerization with multisite catalysts: parameter estimability and experimental design // Polymer Reaction Engineer. 2003. Vol. 11. №3. P. 563-588.

4. Walter E., Pronzato L. Identification of Parametric Models from Experimental Data. Communications and Control Engineering Series. London, UK: Springer, 1997. 413 pp.

5. Brun R., Reichert P., Kunsch H. R. Practical identifiability analysis of large environmental simulation models // Water Resour Res. 2001. Vol. 37. №.4. P. 1015-1030.

6. Быков В. И. Прямые и обратные задачи в химической кинетике. М.: Наука, 1993. 288 c.

7. Губайдуллин И. М. Информационно-аналитическая система решения многопараметрических обратных задач химической кинетики: дис. ... д-ра физ.-мат. наук. Уфа, 2012. 243 с.

8. Turanyi T. Sensitivity analysis of complex kinetic systems. Tools and applications // J Math Chem. 1990. Vol. 5. №3. P. 203-248.

9. Saltelli A., Tarantola S., Campolongo F. Sensitivity analysis as an ingredient of modeling // Statist. Sci. 2000. Vol. 15. №4. P. 377-395.

10. Нурисламова Л. Ф., Губайдуллин И. М. Численный анализ идентифицируемости параметров математической модели химической реакции // Вычислительные методы и программирование. 2018. T. 19. №3. C. 282-292.

11. Nurislamova L. F., Gubaydullin I. M. Mechanism reduction of chemical reaction based on sensitivity analysis: development and testing of some new procedure // Journal of Mathematical Chemistry. 2017. Vol. 55. №9. P. 1779-1792.

12. Полак Л. С., Гольденберг М. Я., Левицкий А. А. Вычислительные методы в химической кинетике. М.: Наука, 1985. 280 c.

13. Vajda H., Rabitz E., Walter and Lecourtier Y. Qualitative and quantitative identifiability analysis of nonlinear chemical kinetic-models // Chem. Eng. Commun. 1989. Vol. 83. P. 191-219.

14. Коробов В. И., Очков В. Ф. Химическая кинетика: введение с Mathcad, Maple, MCS. М.: Горячая линия-Телеком, 2019. 384 c.

15. Box GEP, Hunter W. G., MacGregor J. F., Erjavec J. Some problems associated with the analysis of multiresponse data // Technometrics. 1973. Vol. 15. P. 33-51.

16. Rodriguez-Fernandez M., Egea J. A., Banga J. R. Novel me-taheuristic for parameter estimation in nonlinear dynamic biological systems // BMC Bioinformatics. 2006. Vol. 7. №.483.

17. DAISY Differential Algebra for Identifiability of Systems. URL: https://daisy.dei.unipd.it. (дата обращения: 04.04.2022).

Поступила в редакцию 04.04.2022 г.

ISSN 1998-4812

BecTHHK EamKHpcKoro yHHBepcHTeTa. 2022. T. 27. №2

305

DOI: 10.33184/bulletin-bsu-2022.2.8

ANALYSIS OF THE IDENTIFICABILITY OF KINETIC MODEL PARAMETERS OF CHEMICAL REACTIONS

© L. F. Safiullina1*, I. M. Gubaydullin2

1Bashkir State University 32 Zaki Validi Street, 450076 Ufa, Republic of Bashkortostan, Russia.

2Institute of Petrochemistry and Catalysis, Ufa Federal Research Center of RAS 131 Oktyabrya Avenue, 450075 Ufa, Republic of Bashkortostan, Russia.

*Email: safiullinalf@gmail.com

Kinetic models of chemical reactions usually consist of ordinary differential equations with many unknown parameters. Some of these parameters are often practically unidentifiable, that is, their values cannot be unambiguously determined from the available data. Therefore, before defining model parameter values, it is important to characterize the subset of identified parameters and their interactions. The ability to successfully estimate model parameters depends on the structure of the model, the number of parameters, and the set of experimental data. A set of parameters which evaluation is impossible may be adjusted by applying nominal values, or the model may be modified in order to exclude those parameters at all. Alternatively, the developer of the model may propose additional experiments that may increase the number of identified parameters. The authors of the article present methods to determine a subset of identifiable model parameters. Practical application of the method is shown on the example of a model reaction of alpha-pinene isomerization. Using the DAISY software package for the analysis of the structural identifiability of the model, the authors made sure that the mathematical model is identifiable. The sensitivity analysis of all five parameters of the mathematical model was carried out by the orthogonal method and the eigenvalue method. According to the results of sensitivity and identifiability analyses, the values of model parameters satisfactory correspond to the values of experimental data.

Keywords: sensitivity analysis, identifiability analysis, kinetic model.

Published in Russian. Do not hesitate to contact us at bulletin_bsu@mail.ru if you need translation of the article.

REFERENCES

1. Miao H., Xia X., Perelson A. S., We H. SIAM Rec. Soc. Ind. Appl. Math. 2011. Vol. 53. No. .1. Pp. 3-39.

2. Saccomani M. P., Thomaseth K. Complexity. 2018. Pp. 1-10.

3. Yao K. Z., Shaw B. M., Kou B., McAuley K. B., Bacon D. W. Polymer Reaction Engineer. 2003. Vol. 11. No. 3. Pp. 563-588.

4. Walter E., Pronzato L. Identification of Parametric Models from Experimental Data. Communications and Control Engineering Series. London, UK: Springer, 1997.

5. Brun R., Reichert P., Kunsch H. R. Water Resour Res. 2001. Vol. 37. No. .4. Pp. 1015-1030.

6. Bykov V. I. Pryamye i obratnye zadachi v khimicheskoi kinetike [Direct and inverse problems in chemical kinetics]. Moscow: Nauka, 1993.

7. Gubaidullin I. M. Informatsionno-analiticheskaya sistema resheniya mnogoparametricheskikh obratnykh zadach khimicheskoi kinetiki: dis. ... d-ra fiz.-mat. nauk. Ufa, 2012.

8. Turânyi T. J Math Chem. 1990. Vol. 5. No. 3. Pp. 203-248.

9. Saltelli A., Tarantola S., Campolongo F. Statist. Sci. 2000. Vol. 15. No. 4. Pp. 377-395.

10. Nurislamova L. F., Gubaidullin I. M. Vychislitel'nye metody i programmirovanie. 2018. Vol. 19. No. 3. Pp. 282-292.

11. Nurislamova L. F., Gubaydullin I. M. Journal of Mathematical Chemistry. 2017. Vol. 55. No. 9. Pp. 1779-1792.

12. Polak L. S., Gol'denberg M. Ya., Levitskii A. A. Vychislitel'nye metody v khimicheskoi kinetike [Computational methods in chemical kinetics]. Moscow: Nauka, 1985.

13. Vajda H., Rabitz E. Chem. Eng. Commun. 1989. Vol. 83. Pp. 191-219.

14. Korobov V. I., Ochkov V. F. Khimicheskaya kinetika: vvedenie s Mathcad, Maple, MCS [Chemical kinetics: Introduction with Mathcad, Maple, MCS]. Moscow: Goryachaya liniya-Telekom, 2019.

15. Box GEP, Hunter W. G., MacGregor J. F., Erjavec J. Technometrics. 1973. Vol. 15. Pp. 33-51.

16. Rodriguez-Fernandez M., Egea J. A., Banga J. R. BMC Bioinformatics. 2006. Vol. 7. No. .483.

17. DAISY Differential Algebra for Identifiability of Systems. URL: https://daisy.dei.unipd.it. (data obrashcheniya: 04.04.2022).

Received 04.04.2022.

i Надоели баннеры? Вы всегда можете отключить рекламу.