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

УДК 576/8

Vladislav A. Kholodnov1, Denis A. Krasnoborodko1, Roman Yu. Kulishenko1, Marina Yu. Lebedeva2

NEW IDENTIFICATION METHOD OF MODELS OF COMPLEX CHEMICAL REACTIONS

St Petersburg State Institute of Technology (Technical University), Moskovskiy Pr., 26, St Petersburg, Russia National Research University "Moscow Power Engineering Institute" in Smolensk Energeticheskiy str. 1, Smolensk, 214013, Russia. e-mail: holodnow@yandex.ru

To study the mechanisms of complex chemical reactions and determine their kinetic parameters, the article proposes a method for reduction of optimization parameters. Using the proposed method, the parametric identification of a complex chemical reaction is considered: a set of kinetic constants is found that describes the experimental curves of the concentration-time dependence within this reaction model in the best way. The practical application of the method is illustrated by solving problems of parametric identification of complex chemical reactions using Mathcad. The proposed algorithm can be easily implemented using both traditional and free software.

Keywords: optimization, Mathcad, reverse problems of complex chemical reactions kinetics, reduction of optimization parameters, software.

DOI 10.36807/1998-9849-2020-55-81-91-96

Холоднов В.А.1, Краснобородько Д.А.1, Кулишенко Р.Ю.1 Лебедева М.Ю.2

НОВЫЙ МЕТОД ИДЕНТИФИКАЦИИ МОДЕЛЕЙ СЛОЖНЫХ v ХИМИЧЕСКИХ РЕАКЦИЙ

Санкт-Петербургский государственный технологический институт (технический университет), Московский пр. 26, Санкт-Петербург, Россия Национальный исследовательский университет «МЭИ» в г. Смоленске, Энергетический проезд, 1, г. Смоленск, 214013, Россия. e-mail: holodnow@yandex.ru

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

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

Дата поступления - 7 августа 2020 года

Введение

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

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

ции, разрабатываемый в ННГУ им. Лобачевского. Метод использует редукцию размерности на основе кривых Пеано и информационно-статистический подход, дополненный схемой построения множественных отображений (вращаемые развёртки), позволяющих эффективно использовать сотни процессоров [3].

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

Постановка задачи

Моделирование кинетики сложных химических реакций состоит из двух взаимосвязанных стадий [4].

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

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

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

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

Прямая кинетическая задача для закрытой системы представляет собой задачу Коши для системы обыкновенных дифференциальных уравнений:

rfc,

...,Сп ,к1,к2,..., кщ),

или в виде

С начальными условиями ( = О, С\(0) = с;0 где С, - концентрация (мольная доля) /-го вещества, участвующего в химической реакции; / = 1,2, ... п; к -константа скорости /ой стадии реакции; j = 1,2, ... т; п - количество веществ; t- время; кр- предэкспонента р ой стадии реакции; р= 1,2, ... т0; р - энергия активации предэкспонента р ой стадии реакции; р = 1,2, ...

то.

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

Моделирование кинетики сложных химических реакций тесно связано с решением прямой и обратной задачи.

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

г-/,, ,, \ _ ^дрдстящг^&ияитг\2

* ■■■ ".V _ ¡-■■■-=- —:■■=. 1 -_____:.■ ,

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

Математическая постановка задачи имеет следующий вид:

найти I: = г: (

с учетом ограничений на параметры оптимизации : .

В качестве параметров оптимизации используются

или к} - константа скорости / - ой стадии реакции; / = 1,2, ...,т или к}0 - предэкспонента /-ой

стадии реакции; / = 1,2, ...,т0 и Е,- - энергия активации / - ой стадии реакции;/ = 1,2,..., тО.

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

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

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

Для решения такого рода задач оптимизации, связанных с решением на каждом шаге оптимизации систем дифференциальных уравнений, нами предлагается следующий оригинальный метод-редукция параметров оптимизации. В объект оптимизации - систему дифференциальных уравнений добавляются дифференциальные уравнения по ] -му параметру оптимизации и- в виде: ¡Ь^/ск = 0. Начальные условия для параметра ц меняются на каждом шаге оптимизации, что

позволяет существенно упростить решение задачи оптимизации.

Приведем использование данного метода для решения традиционной задачи в химической технологии [5].

Определение кинетических констант для системы химических реакций при одной температуре

Дана система химических реакций:

ki

+B -

->C

Зависимость концентраций веществ от времени можно представить в виде системы дифференциальных уравнений:

¿а ,

— = -к,а

а

— = ка - кЪ ¿1 1 2

¿С = к2Ъ

А 2

В системе кинетических уравнений малыми буквами обозначены концентрации соответствующих веществ. Известны значения концентраций веществ в начальный момент времени:

при I = 0 а0 = 100 Ь0 = 0, с0 = 0.

По этим данным требуется определить параметры кЛ и к2.

С помощью вычислительного эксперимента при заданных значениях к1 и к2, решая систему рассмотренных выше дифференциальных уравнений в рамках МаШСАЭ, получаем «точные» значения кон-

центраций по всем компонентам, в моменты времени t [0,10].

Задание начального приближения неизвестных констант:

к1:= 0.5 k2:= 0.5 п := 0 tk := 10 п := 100

Задание начальных концентраций:

А 00 ^ 0 0 k1 к 2 У

Точные значения концентраций А, В, С в моменты времени 2, 4, 6, 8, 10:

c exp t :=

Последние два значения вектора x представляют параметры оптимизации k1 и k2 .

С помощью встроенной функции rnd, из «точных данных» получаем «экспериментальные данные». Используя найденные значения «экспериментальных данных», с помощью встроенной функции Minimize, решаем задачу оптимизации по нахождению значений предэкспонент и энергий активаций.

«Экспериментальные» значения концентраций по результатам вычисленного эксперимента: t ;= 10 i := 1..5 j := 1..3

cexp^. := cexpt ^ + :j~'cexPttJ • (1 - rnd(1)) • sign(0.5 - rnd(1))

' 30.3 49.64 20.35N

9.0l 42.l2 48.88

2.8 37.8 69.73

0.9 l6.46 82.73

v 0.26 9.l 90.7 ,

'33.326 51.7 20.7llN

9.27l 45.947 48.936

2.522 28.907 73.56

0.88 l6.256 86.54l

ч 0.256 9.097 84.045 y

c exp =

Программа интегрирования системы дифференциальных уравнений:

X4 • X1

D(t, x) :=

г := Rkadapt(x, п, tk, п, В)

И := г« А := г<2> В := 3 С := г<4 12 := 2 • . А ехр := с ех^п В ехр := с ехр; 2

С ехр,. := с ехР,,з

Последние два значения вектора D(t,x) характеризуют равенство 0 соответствующих производных.

На рис. 1 представлено сравнение рассчитанных значений с экспериментальными.

Аехр

Вехр 60

Сехр

/

! ~-—i

t2.t1t2.t1.t2. И

Рис. 1. Экспериментальные и рассчитанные значения концентраций А,В, С

тодом:

Программа решения задачи предлагаемым ме-

TOL := 0.001 CTOL := 0.001 x0 ^ x x04 ^ k1 x05 ^ k 2 i ^ 1 S ^ 0

for t e 2,4,6,8,10

Z ^ AdamsBDF(x0,0, t, n, D)

OPT (kl, k2) :=

Ql ^ (Z 4 - c expu)2

kl ^ (kl

I := Minimize(OPT, kl, k2) I k 2 I I k 2

Q2 ^ (Zn,3 - c exp,2)2 Q3 ^ (Zn,2 - c expij )2 S ^ S + Q1 + Q2 + Q3 i ^ i +1 Q ^ s

OPT(k1, k2) = 1.996 x103 Программа решения задачи оптимизации:

Given

r 0.611 ^ , v0.299) OPT (k1, k2) = 100.062 Результат решения задачи: x0 := x

Z1 := AdamsBDF (x0,0,tk, n, D) t1:= Z1« A := Z1<2> B := Z1<3> C := Z1(4> 12i := 2 • i A exp := c exp; j

B exp;. := c expi,2 Cexp := cexp, 3

На рис. 2 представлено сравнение экспериментальных и рассчитанных значений при найденных оптимальных значениях кЛ и к2.

^Сл * *

Х5 ' С2

V

/

C + D

— D + E

2 • T

2 • T

С помощью вычислительного эксперимента решая систему рассмотренных выше дифференциальных уравнений в рамках Ма^САЭ, получаем «точные» значения концентраций по всем компонентам, в том числе и по ключевым компонентам А, С, Е при различных температурах 283К, 293К, 303К, 313К в моменты времени I [0,20].

С помощью встроенной функции rnd, из «точных данных» получаем «экспериментальные данные». Используя найденные значения «экспериментальных данных» при 4-х температурах 283К, 293К, 303К и 313К, с помощью встроенной функции Minimize, решаем задачу оптимизации по нахождению значений предэкспонент и энергий активаций. Определяем такие значения предэкспонент и энергий активаций, чтобы сумма квадратов относительных отклонений «экспериментальных» значений концентраций при 4-х температурах в 5-и моментах времени 4, 8, 12, 16, 20 по 3 ключевым компонентам A, C и E от расчетных была минимальной. Целевая функция имеет вид: F(fe10,ir20,fe30,E1„E2„£3)

Рис. 2. Экспериментальные и рассчитанные значения концентраций А,В, С при оптимальных к1 и к2

Определение кинетических констант для системы химических реакций при разных температурах

Суммарная реакция преобразования одновалентного таллия Т1+ (элемента третьей группы таблицы Д.И. Менделеева) в трехвалентный Т13+ при взаимодействии с трехвалентным кобальтом Со3+ имеет промежуточный продукт (двухвалентный таллий Т?+) и происходит в две стадии, причем первая из них обратимая:

а+в к1 > с+Б

£пТ-=1 Sr=l ¿jfc

' ? , (—-— Т.-пщгу тш

1 ■ Пчкггг.

а+в < в+с -

В системе приняты обозначения: А~ Т+, В- СО+, С- Т2+, й~ Со2+, Е-Т3+.

Кинетика реакции Т++2 СО+ — Т++2СО+ в периодическом реакторе описывается системой из пяти дифференциальных уравнений:

— = ке ■ ё - к,а ■ Ъ ¿г 2 1

— = ке ■ ё - к,а ■ Ъ - кЪ ■ е

¿г 2 1 3

— = к,а ■ Ъ - ке ■ ё - к,Ъ ■ е

¿г 1 2 3

-= к, а ■ Ъ - ке ■ ё + кф ■ е

¿г 1 2 3

аЕ ьи

— = к,Ъ ■ е

йг 3

При моделировании были приняты следующие значения концентраций веществ в начальный момент времени: при 0 = 0, а0 = Ь0 = 1, с0 = с<0 = е0 =0.

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

1ЛГТЛ П. 14000 ч

к1(Т ) = ехр(23 ———), 2 • Т

... 16000, 18000 к 2(Т) = ехр(24--), к 3(Т) = ехр(30--)

Результат решения этой задачи представлен на рис. 3-6.

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

В качестве исходной точки поиска была взята

точка:

1п к10 = 22, 1п к10 = 25, 1п к10 = 28,

Е1 = 11100. Е2 = 13000, Е1 = 17000 со значением целевой функции Е = 0.628.

Задание начального приближения неизвестных констант:

х1:= 23 х2 := 14000 х3 := 24 х4:= 16000 х5:= 30 х6:= 18000

tn := 0 tk := 20 h := 4

100 • (tk - tn) h

Точные значения кон центраций А, В, С в моменты времени 2, 4, 6, 8, 10 при температурах 283К, 293К, 303К, 313К:

(0.393 0.257 0.35 ^

T283 :=

0.433 0.161 0.451 0.117

0.407 0.432 0.461 0.092 0.447

0.468 0.076 0.456

(0.479 0.052 0.469^

0.489 0.027 0.484

T303:= 0.493 0.018 0.489

0.495 0.014 0.492

0.496 0.011 0.493

(0.45 0.122 0.428^

0.472 0.068 0.46

T293:= 0.481 0.047 0.473

0.485 0.036 0.479

0.488 0.029 0.483^

(0.492 0.022 0.487^

0.496 0.011 0.493

T313:= 0.497 0.0073 0.496

0.498 0.0055 0.497

0.498 0.0044 0.497

ct := stack(T283,T293, T303, T13) «Экспериментальные» значения концентраций по результатам вычисленного эксперимента: t := 2 i := 1..20 j := 1..3

c exp,.j. := ctij + щ • ctij • (1 - rnd (1)) • sign(0.5 - rnd (1))

Программа интегрирования системы дифференциальных уравнений на отрезке tn, tk и вычисления суммы квадратов относительных отклонений концентраций:

( c c \ - exp(c6 - • c1 • c2 + exp(c8 - уМ • c3 • c4

2 • ciO 2 • ciO

- exp(c6 - Y^-) * cl * c2 + eXP(C8 - * C3 * C4 - exp(cl0 — ^^ * C2 ' C3

'(Сю — ^ 2 * c,.

exp(C6 - * Cl * C2 - exp(C8 - -^-l 2 * Cl2 2 * Cl2

*c *c

* C * C

exp(c6 - y-M * Cl * C2 - exp(c8 - yM * C3 * C4 + exp(cl0 - Y^-') * C2 * C3

D(t' C) := exp(cl0 - Т^Ч * C2 * C3 2 * Cl2

0 0 0 0 0 0 . 0

Обращение к программе оптимизации:

TOL := 0.001 CTOL := 0.001 j ^ 0 k ^ 0 S ^ 0

100 • (tk - tn)

n <-----

h

forT e 283,293,303,313

f1 ^ 1 0 0 0 x1 x2 x3 x4 x5 x6 T __

Z ^ rkfixed(c, tn, tk, n, D) for i e 1..5

opt(xl, x2, x3, x4, x5, x6):

^ Zi*l00 +l,2 C expi+j,ly

QC ^ (

c expi+j,l Zi*l00 +l,4 — C expi+j,2

c exR+

"c exR+

QE ^_^ — i lUU +l,o - —ri+j,3

c exR+j,3 j ^ j +1

opt (x1, x2, x3, x4, x5, x6) = 616.431 Программа решения задачи оптимизации:

Given

x1 > 22 x1 < 25.6 x2 > 12000 x2 <15500 x3 > 23 x3 < 27 x4 > 13500 x4 < 17000 x5 > 29 . 1 x5 < 33 x6 > 16000 x6 < 19000

(xl > ( x l > '24.2494

x2 x2 l3414

x3 x3 25.005

x4 := Minimize(opt, xl, x2, x3, x4, x5, x6) =

x4 l492l

x5 x5 30.75l

vx6 ) ч x6 > l7286

c0 :=

opt (x1, x2, x3, x4, x5, x6) = 0.101 Результат решения задачи при температуре 283К: T := 283

n ^ 1 0 0 0 x1 x2 x3 x4 x5 x6

z := Rkadapt(c0, tn, tk, n, D)

i := 1..5 t1:= Z1 A := Z3 C := Z4) E := Z6 12,:= 4• i Aexp:= T283,.j

C exp := T283,-2 Eexp := T2 83,3

На рис. 3-6 представлены экспериментальные и рассчитанные значения концентраций компонентов А, С, Е при различных температурах.

Аехр • • •

А

Сехр 0.1 ■ ■ ■

С

Еехр 0.

А А

Е

4 у ч

''i i— ---------,

8 12 t2.t1.t2.t1,12,«

Рис. 3. Экспериментальные и рассчитанные значения при температуре 283К

S

Рис. 4. Экспериментальные и рассчитанные значения при температуре 293К

Рис. 5. Экспериментальные и рассчитанные значения при температуре 303К

Рис. 6. Экспериментальные и рассчитанные значения при температуре 313К

Выводы

Предложенный метод позволяет успешно решать задачу идентификации сложных химических реакций.

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

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

Литература

1. Л.С. Полак. Применение вычислительной математики в химической и физической кинетике, М: Наука, 1969, 282 с.

2. Гольдштейн АЛ.Оптимизация в среде MATLAB: учеб. пособие Пермь: Изд-во Перм. Нац. Ис-след. Политехн. Ун-та, 2015. 192 с.

3. Тихонова М.В., Рябов В.В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислит. методы и программирование. 2011. Т 12. Вып. 1. С. 137-145

4. Холоднов В.А., Боровинская ЕС, Андреева В.П., Черемисин В.И. Решение задач нелинейного программирования на основе градиентных методов с использованием системы компьютерной математики MathCAD: метод. указания СПб.: СПбГТИ (ТУ), 2010. 69 с.

5. Холоднов В .А., Дьяконов В.П., Фонарь В.В., Кулишенко Р.Ю, Ананченко И.В. Системный анализ и принятия решений. Технология вычислений в системе компьютерной математики Mathcad: учеб. пособие. СПб.: СПбГТИ (ТУ), 2012. 153 с.

References

1. LS. Polak. Primenenie vychislitel'noj ma-tematiki v himicheskoj i fizicheskoj kinetike, M: Nauka, 1969, 282 s.

2. Gol'dshtejn A.L. Optimizaciya v srede MATLAB: ucheb. posobie Perm': Izd-vo Perm. Nac. Issled. Politekhn. Un-ta, 2015. 192 s.

3. Tihonova M.V, Ryabov V.V. Primenenie in-deksnogo metoda global'noj optimizacii pri reshenii obratnyh zadach himicheskoj kinetiki // Vychislit. metody i programmirovanie. 2011. T 12. Vyp. 1. S. 137-145

4. Holodnov V.A., BBorovinskaya ES, Andreeva V.P, Cheremisin V.I. Reshenie zadach nelinejnogo pro-grammirovaniya na osnove gradientnyh metodov s ispol'zovaniem sistemy komp'yuternoj matematiki MathCAD: metod. ukazaniya SPb.: SPbGTI (TU), 2010. 69 s.

5. Holodnov V.A., D'yakonov V.P, Fonar' V.V, Kuiishenko R.YU, Ananchenko I.V. Sistemnyj analiz i prinyatiya reshenij. Tekhnologiya vychislenij v sisteme komp'yuternoj matematiki Mathcad: ucheb. posobie. SPb.: SPbGTI (TU), 2012. 153 s.

Сведения об авторах

Холоднов Владислав Алексеевич, д-р техн. наук, профессор, каф. системного анализа; Vladislav A. Kholodnov, Dr Sci. (Eng.), Professor, Department of Systems Analysis holodnow@yandex.ru

Краснобородько Денис Александрович, канд. техн. наук, доцент, каф. системного анализа; Denis A. Krasnoborodko, Ph.D. (Eng.), Associate Professor, Department of Systems Analysis

Кулишенко Роман Юрьевич, канд. техн. наук, доцент, каф. системного анализа; Roman Yu. Kulishenko Ph.D. (Eng.), Associate Professor, Department of Systems Analysis

Лебедева Марина Юрьевна, канд. техн. наук, доцент каф. информационных технологий в экономике и управлении; Marina Yu. Lebedeva, Ph.D. (Eng.), Assistant Professor, Department of management and information technology in the economy, marilieb@yandex.ru