Интервальный подход при решении задач кинетики простых химических реакций Текст научной статьи по специальности «Математика»

Вычислительные технологии

Том 2, № 1, 1997

ИНТЕРВАЛЬНЫЙ ПОДХОД ПРИ РЕШЕНИИ ЗАДАЧ КИНЕТИКИ ПРОСТЫХ ХИМИЧЕСКИХ РЕАКЦИЙ*

В. М. Белов, В. А. Суханов, Е. В. Лагуткина Алтайский государственный технический университет

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

Обзор литературы показывает, что в основу интервального оценивания параметров эмпирических зависимостей рядом авторов положена работа Л. В. Канторовича [1], в которой проведена постановка задачи определения возможных границ параметров Хк известной функциональной зависимости в форме

п

Уг - $ < ^ С^к Хк < Уг + « = 1,Ш.

к= 1

Решение этой задачи автор [1] предложил искать, используя методы линейного программирования.

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

До настоящего времени среди интервальных методов оценивания параметров физико-химических моделей серьезное развитие получил лишь метод [3-14], сочетающий принцип выравнивания по П. Л. Чебышеву с методологией Л. В. Канторовича и успешно примененный в задачах обработки данных кинетических измерений (благодаря указанным публикациям в литературе появился термин "нестатистический подход").

При математической интерпретации геохимических данных близкие подходы развивал И. К. Карпов, подробно рассмотревший задачу выпуклого программирования с ограничениями в исходной геологической информации [15-19]. Г. А. Коковин и В. А. Титов [2, 20-22] показали перспективность применения алгоритмов нестатистического (интервального) оценивания к задачам поиска констант равновесия в ситуациях, когда в качестве погрешностей взяты предельно возможные ошибки измерений.

Первые работы по приложению нестатистических интервальных методов (будем их так называть) в аналитической химии опубликовали Г. И. Рузайкин [23] и авторы настоящей статьи [24-29].

*© В.М.Белов, В.А.Суханов, Е. В. Лагуткина, 1997

Среди работ, имеющих отношение к нестатистическому оцениванию параметров эмпирических зависимостей, дополнительно выделим работы А. П. Вощинина [30, 31].

1. Определение констант скорости необратимых реакций первого порядка

В химической кинетике обобщенное стехиометрическое уравнение необратимой реакции первого порядка обычно записывают в виде [32, 33]

т

аА — Ь Б, (1)

3=1

где А — исходное вещество, Б3- — продукты химической реакции, а и — стехиомет-рические коэффициенты. Зависимость концентрации вещества А от времени определяют решением задачи Коши для обыкновенных дифференциальных уравнений типа [32, 33]

адм = -К *[А], [А] (¿о) = [А]о, (2)

где К * — константа скорости необратимой реакции первого порядка, [А] — концентрация вещества А.

1.1. Постановка интервальной задачи

Константу скорости необратимой химической реакции, согласно задаче (2), будем искать по формуле [32, 33]

К * = (1п([А]о/[А]г))/(^ - ¿о), (3)

где [А]о, [А]г — соответственно концентрации вещества А при £ = ¿о и более позднем моменте времени Через А*(г) и Т*(г) обозначим истинные значения концентрации вещества А и времени £ при г-м измерении. Для нашей задачи интервальная модель обработки измерений задается следующими соотношениями:

Т(г) - г! < Т*(г) < Т(г) + гь (4)

А(г) - £2 < А*(г) < А(г) + £2, (5)

А(0) = А*(0), Т(0) = Т*(0), г е Т,Ж (6)

В формулах (4)-(6) через А(г) и Т(г), г е Т,Ж обозначены экспериментальные значения концентрации [А] и времени £ при г-м измерении. Общее число проведенных измерений обозначено через N.

1.2. Решение интервальной задачи

Используя формулу (3) и введенные выше обозначения, получим

К * = (1п(А*(0)/А*(г)))/(Т *(г) - Т*(0)), г е Т,^

(7)

Из формул (4)-(7) имеем

1п(А(0)/(А(г) + е2)) < К* < 1п(А(0)/(А(г) - ¿2)) Т(г) - Т(0)+ ¿1 - - Т(г) - Т(0) - ¿1 " и

Введем обозначения:

К + (г) = (1п(А(0)/(А(г) - ¿2)))/(Т(г) - Т(0) - ¿1), (9)

К-(г) = (1п(А(0)/(А(г) + ¿2)))/(Т(г) - Т(0) + ¿1), (10)

К- = шах К-(г), К + = тт К+(г), (11)

Тогда выражение (8) можно переписать в виде

К-(г) < К < К + (г), г е 17Ж (12)

Из соотношений (11), (12) следует, что должно выполняться неравенство

К- < К* < К+ ^ К* е [К-,К+]. (13)

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

Множество вещественных чисел отрезка [К-, К+], следуя работе [33], назовем множеством неопределенности неизвестной константы скорости К*. Центр тяжести этого множества примем в качестве оценки неизвестного нам значения константы скорости реакции. Таким образом,

К = (К- + К+)/2. (14)

Из (13), (14) имеем

|К - К*| < е для £ = 0.5 ■ (К + - К-). (15)

Число е оценивает сверху погрешность, с которой определяем константу скорости односторонней химической реакции К *. При обработке результатов реального эксперимента возможна ситуация, когда получится соотношение К- > К +. Это означает, что неверно оценены значения ошибок е1 и е2 измерения величин Т*(г) и А*(г).

1.3. Алгоритмическое обеспечение

Для алгоритмизации расчетной методики примем следующие обозначения:

5 (г) = К (1) + К (2) + ... + К (г - 1) + К (г), (16)

где К (г) = [1п(А(0)/А(г))]/[Т (г) - Т (0)]);

P+ (i) = min K+(i), P- = max K-(i) i £ 1, N. (17)

1<j<i 1<j<i

Из (16), (17) получаем соотношения

P+(0) = 0, P-(0) = 0, S (0) = 0, (18)

S(i) = S(i - 1) + K(i),

P+(i) = min{P+(i - 1),K+(i)}, (19)

P-(i) = max{P-(i - 1),K-(i)}, i £ 17N.

Kcp = S (N )/N, K + = P+(N), K- = P -(N), K = (K+ + K-)/2. (20)

С использованием формул (16)-(20) легко построить пошаговый алгоритм решения интервальной задачи.

1.4. Пример вычислений

В качестве примера вычислений была взята реакция распада гексафенилэтана на два свободных радикала [33, 34]:

(Об^зО - 0(06^)3 2(06^)03 ■. (21)

В табл. 1 представлена исходная информация по кинетике реакции (21).

Таблица 1. Разложение гексафенилэтана при 0°С в смеси 95% толуола и 5% анилина

г, мин 0 0.50 1.05 2.20 3.65 5.5 7.85 9.45 14.75

Концентрация [А] гексафенилэтана, моль/л 0.1000 0.0934 0.0867 0.0733 0.0600 0.0465 0.0334 0.0265 0.0134

г 1 2 3 4 5 6 7 8 9

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

Т-(г) = Т(г) - еь Т+(г) = Т(г) + гь (22)

А-(г) = А(г) - £2, А+(г) = А(г) + £2, (23)

Т-(0) = Т+(0) = Т (0), А-(0) = А+(0) = А(0), (24)

К (г) = [1п(А(0)/А(г))]/[Т (г) - Т (0)], г € (25)

Из табл. 2 видно, что К- > К+ при е1 = 0.005 и е2 = 0.00005, а также при е1 = 0.005 и е2 = 0.0001. Следовательно, принятые в этих случаях предположения относительно ошибок определения времени и концентрации не соответствуют реальности. Далее в приведенных в табл. 2 данных такого противоречия не обнаружено.

2. Одновременное оценивание константы скорости необратимой реакции первого порядка и начальной концентрации реагента

2.1. Постановка интервальной задачи

Рассмотрим ту же реакцию распада гексафенилэтана, что и в разделе 1, со стехиомет-рическим уравнением, задаваемым выражением (1). Дифференциальное уравнение, отражающее связь концентрации вещества А и времени протекания реакции, соответствует задаче (2).

Решение (2) относительно концентрации вещества А имеет вид

[А](*) = [А]о ехр(-К(г - ¿с)), г > ¿с.

(26)

Таблица 2. Оценка константы скорости реакции разложения гексафенилэтана

Время Концентрация Константы скорости реакции

г Т" (г) Т (г) Т +(г) А-(г) А(г) А+ (г) К" (г) К (г) К+(г)

е1 = 0.005, £2 = 0.00005

0 0.0 0.0 0.0 0.1 0.1 0.1 - - -

1 0.495 0.50 0.505 0.09335 0.0934 0.09345 0.13415 0.13656 0.13902

2 1.045 1.05 1.055 0.08665 0.0867 0.08675 0.13473 0.13592 0.13712

3 2.195 2.20 2.205 0.07325 0.0733 0.07335 0.14056 0.14119 0.14182

4 3.645 3.65 3.655 0.05995 0.0600 0.06005 0.13953 0.13995 0.14037

5 5.495 5.50 5.505 0.04645 0.0465 0.04655 0.13890 0.13922 0.13954

6 7.845 7.85 7.855 0.03335 0.0334 0.03345 0.13942 0.13970 0.13998

7 9.445 9.45 9.455 0.02645 0.0265 0.02655 0.14026 0.14053 0.14081

8 14.745 14.75 14.755 0.01335 0.0134 0.01345 0.13597 К- = 0.14056 0.13627 Кср = 0.13867 К = 0.13856 0.13657 К+ = 0.13657

£1 = 0.005, £2 = 0.0001

0 0.0 0.0 0.0 0.1 0.1 0.1 - - -

1 0.495 0.50 0.505 0.0933 0.0934 0.0935 0.1331 0.1366 0.1401

2 1.045 1.05 1.055 0.0866 0.0867 0.0868 0.1342 0.1359 0.1377

3 2.195 2.20 2.205 0.0732 0.0733 0.0734 0.1403 0.1412 0.1421

4 3.645 3.65 3.655 0.0599 0.0600 0.0601 0.1393 0.1400 0.1406

5 5.495 5.50 5.505 0.0464 0.0465 0.0466 0.1387 0.1392 0.1397

6 7.845 7.85 7.855 0.0333 0.0334 0.0335 0.1392 0.1397 0.1402

7 9.445 9.45 9.455 0.0264 0.0265 0.0266 0.1401 0.1405 0.1410

8 14.745 14.75 14.755 0.0133 0.0134 0.0135 0.1357 К- = 0.1403 0.1363 Кср = 0.1387 К = 0.1385 0.1368 К+ = 0.1368

£1 = 0.005, £2 = 0.0005

0 0.0 0.0 0.0 0.1 0.1 0.1 - - -

1 0.495 0.50 0.505 0.0929 0.0934 0.0939 0.1246 0.1366 0.1488

2 1.045 1.05 1.055 0.0862 0.0867 0.0872 0.1298 0.1359 0.1421

3 2.195 2.20 2.205 0.0728 0.0733 0.0738 0.1378 0.1412 0.1446

4 3.645 3.65 3.655 0.0595 0.0600 0.0605 0.1375 0.1400 0.1424

5 5.495 5.50 5.505 0.0460 0.0465 0.0470 0.1372 0.1392 0.1413

6 7.845 7.85 7.855 0.0329 0.0334 0.0339 0.1377 0.1397 0.1417

7 9.445 9.45 9.455 0.0260 0.0265 0.0270 0.1385 0.1405 0.1426

8 14.745 14.75 £1 = 0 14.755 £2 0.0129 = 0.0005 0.0134 0.0139 0.1337 К- = 0.1385 0.1363 Кср = 0.1387 К = 0.1387 0.1389 К+ = 0.1389

0 0.0 0.0 0.0 0.1 0.1 0.1 - - -

1 0.50 0.50 0.50 0.0929 0.0934 0.0939 0.1259 0.1366 0.1473

2 1.05 1.05 1.05 0.0862 0.0867 0.0872 0.1304 0.1359 0.1414

3 2.20 2.20 2.20 0.0728 0.0733 0.0738 0.1381 0.1412 0.1443

4 3.65 3.65 3.65 0.0595 0.0600 0.0605 0.1377 0.1400 0.1422

5 5.50 5.50 5.50 0.0460 0.0465 0.0470 0.1373 0.1392 0.1412

6 7.85 7.85 7.85 0.0329 0.0334 0.0339 0.1378 0.1397 0.1416

7 9.45 9.45 9.45 0.0260 0.0265 0.0270 0.1386 0.1405 0.1426

8 14.75 14.75 14.75 0.0129 0.0134 0.0139 0.1338 К- = 0.1386 0.1363 Кср = 0.1387 К = 0.1387 0.1388 К+ = 0.1388

£1 = 0, £2 = 0.001

0 0.0 0.0 0.0 0.1 0.1 0.1 - - -

1 0.50 0.50 0.50 0.0924 0.0934 0.0944 0.1153 0.1366 0.1581

2 1.05 1.05 1.05 0.0857 0.0867 0.0877 0.1250 0.1359 0.1470

3 2.20 2.20 2.20 0.0723 0.0733 0.0743 0.1350 0.1412 0.1474

4 3.65 3.65 3.65 0.0590 0.0600 0.0610 0.1354 0.1400 0.1446

Нам известны концентрации вещества А с точностью е2 в измеряемые с точностью е1 моменты времени Т*(г). Начало временного отсчета Т*(0) предполагается известным абсолютно точно (Т*(0) = Т(0)). Цель обработки имеющихся массивов (чисел А(г) и Т(г), г е 1, N) эмпирической информации — получение интервальной оценки константы скорости К (К_ < К < К+) и начальной концентрации реагента А (А0- < [А]0 < А0+).

2.2. Решение интервальной задачи

Из формулы (26) следует справедливость соотношения

[А](т2)/[А](т1) = ехр(-К(т2 - Т1)), Т2 > п > ¿о, (27)

откуда можно получить

К = [1п([А](г2)/[А](п))]/(г2 - Т1), Т2 > Т1 > ¿0. (28)

Преобразуя (28), получаем неравенство

1п((А(,?) - £2)/(А(г)+ ¿2)) < К < 1п((А(?)+ ¿2)/(А(г) - ¿2))

Т(?) - Т(г) + 2£1 < < Т(?) - Т(г) - 2£1 ' ( )

где i G 1, N — 1, j > i + 1. Положим

K + (i, j) = [ln((A(j) + e2)/(A(i) — e2))]/(T (j) — T (i) — 2ex), (30)

K -(i, j) = [ln((A(j) — £2)/(A(i) + £2))]/(T (j) — T (i) + 2ei), (31)

где i G 1, N — 1, j > i + 1. Соотношение (29) с учетом (30) и (31) принимает вид

K-(i,j) < K < K +(i,j), i G 1, N — 1, j G i + 1, N. (32)

Введем обозначения

K+ = min min K+(i, j), K- = max max K (i, 7). (33)

1<i<N 1+i<j<N 1<i<N 1+i<j<N

Из формул (32), (33) вытекает неравенство

K- < K < K+. Используя соотношение (26), получаем формулу

[A]o = [A](t)/exp(—K (t — to)), (34)

поэтому

[A](t)/exp(—K-(t — to)) < [A]o < [A](t)/exp(—K+ (t — to)).

Следовательно,

A(i) — e2 r 41 A(i) + e2 —— 2 < [A]o <-trs+w-i—^-TT , i G 1, N.

exp(K-(T(i) — T(0) — ei)) " L Jo" exp(K +(T(i) — T(0) + £i)):

Введем обозначения

A0-(i) = (A(i) — £2)/exp(K-(T(i) — T(0) — ei)), i G 1, N,

(35)

(36)

A0+(i) = (A(i) + £2)/exp(K+(T(i) — T(0) + ei)), i G 1,N,

с учетом которых получим систему неравенств

A0-(i) < [A]o < A0+(i), i G 1, N.

(37)

Обозначим через A0+ и A0 величины:

A0+ = min A0+(i), A0

1<i<N

max A0 (i).

1<i<N

(38)

Из соотношений (37), (38) вытекают неравенства

А0- < [А]0 < А0+.

(39)

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

Список литературы

[1] Канторович Л. В. О некоторых новых подходах к вычислительным методам и обработке наблюдений. Сиб. мат. журн. 3, №5, 1962, 701-709.

[2] Коковин Г. А. и др. Некоторые методологические вопросы математической обработки экспериментальных данных по исследованию равновесий. В "Матем. в химической термодинамике", Новосибирск, 1980, 50-58.

[3] СпивАк С. И., Шмелев А. С. Методологические аспекты определения физико-химических параметров по экспериментальным данным. Там же, 84-91.

[4] СпивАк С. И. и др. Планирование методами математического программирования экспериментов для уточнения параметров кинетических моделей. Кинетика-2. Тез. докл. II Всесоюз. конф. по кинетике катал. реакций, 3, ИК СО АН СССР, Новосибирск, 1975, 36-42.

[5] СпивАк С. И. и др. Методы построения кинетических стационарных реакций. Химическая промышленность №3, 1979, 33-36.

[6] СпивАк С. И. О неединственном решении задач восстановления констант химической кинетики и констант химических равновесий. В "Матем. проблемы химической термодинамики", Наука, Новосибирск, 1980, 63-72.

[7] СпивАк С. И. Информативность эксперимента и проблема неединственности решения обратных задач химической кинетики: Автореф. дис. ... д-ра физ.-мат. наук. Черноголовка, 1984.

[8] Слинько М. Г. О критериях определения параметров кинетических моделей. Кинетика и катализ 13, №6, 1972, 1570-1578.

[9] СпивАк С. И., ТимошЕнко В. И. Применение метода выравнивания по П. Л. Че-бышеву при построении кинетической модели сложной химической реакции. Докл. АН СССР 192, №3, 1970, 580-582.

[10] АсАдуллин Р. М., СпивАк С. И. О критериях определения констант фазовых равновесий. Журн. физ. химии 54, №4, 1980, 890-893.

[11] Ицкович И. А., СпивАк С. И. Анализ применения методов линейного программирования при построении кинетической модели сложной химической реакции. В "Управляемые системы", вып. 4-5, Наука, Новосибирск, 1970, 142-147.

[12] СпивАк С. И. и др. Оценка значимости влияния измерения на кинетическую модель химической реакции. В "Матем. проблемы химии, часть 2", ВЦ СО АН СССР, Новосибирск, 1973, 3-9.

[13] Болдырев В. И., СпивАк С. И. Чебышевские приближения для кинетической модели с дробно-линейной зависимостью от параметров. Там же, 58-65.

[14] АндрушкЕвич М. М. и др. Кинетика процесса окислительного дегидрирования н-бутиленов на хромкальцийникельфосфатном катализаторе. Кинетика и катализ 11, №6, 1970, 1419-1425.

[15] Карпов И. К. и др. Расчет на ЭВМ необратимой эволюции геохимических систем методами оптимального программирования. Геохимия №4, 1973, 603-611.

[16] Карпов И. К. и др. Термодинамика природных мультисистем с ограничивающими условиями. Наука, Новосибирск, 1976.

[17] Карпов И. К. и др. Моделирование природного минералообразования на ЭВМ. Недра, М., 1976.

[18] Карпов И. К. Теоретические основы физико-химического моделирования природных процессов минералообразования на ЭВМ: Автореф. дис. ... докт. геол.-минер. наук. Иркутск, 1978.

[19] Карпов И. К. Физико-химическое моделирование на ЭВМ в геохимии. Наука, Новосибирск, 1981.

[20] Титов В. А., Коковин Г. А. К вопросу о выборе целевой функции при обработке экспериментальных данных по давлению насыщенного пара. В "Матем. проблемы химии, часть 2", ВЦ СО АН СССР, Новосибирск, 1975, 25-34.

[21] Титов В. А. Анализ некоторых физико-химических равновесий с помощью статического метода измерения давления пара. Автореф. дис. ... канд. хим. наук. Новосибирск, 1980.

[22] ТитовА е. Ф., Коковин Г. А. О применении нестатистического подхода к нелинейным обратным задачам. В "Матем. методы химической термодинамики", Наука, Новосибирск, 1982, 170-177.

[23] РузАйкин Г. И. Метод построения по экспериментальным данным сглаженной кривой с контролем за выбором для нее погрешности. В "Матем. методы и ЭВМ в аналит. химии: Тез. докл. Всесоюз. конф.", Наука, М., 1986, 54.

[24] Белов В.М., Суханов В. А., Унгер Ф.Г. Метод центра неопределенности при расчете линейных градуировочных графиков и метрологических характеристик результатов химического анализа. Препр. Ин-та химии нефти, ТНЦ СО АН СССР, №59, Томск, 1989.

[25] Белов В. М. и др. Оценивание параметров линейных физико-химических зависимостей прямоугольником метода центра неопределенности. Изв.вузов. Физика №8, 1991, 35-45.

[26] Белов В. М. и др. Методическое и метрологическое обеспечение рентгенофазового анализа нефтяных асфальтенов на степень кристалличности. Изв. вузов. Химия и хим. технология 36, №5, 1993, 46-50.

[27] Белов В. М. и др. Статистический и нестатистический анализ градуировочных зависимостей в инверсионной вольтамперометрии тяжелых металлов. Там же, №11, 1993, 47-50.

[28] Белов В. М., Суханов В. А., Рыбальченко Т. А. Оценивание параметров линейных физико-химических зависимостей эллипсом в методе центра неопределенности. Там же, 38, №6, 1995, 42-46.

[29] Белов В.М., Суханов В. А., Унгер Ф.Г. Теоретические и прикладные аспекты метода центра неопределенности. Наука, Новосибирск, 1995.

[30] Вощинин А. П., Сотиров Г. Р. Оптимизация в условиях неопределенности. Изд. МЭИ (СССР), Техника (НБР), М.-София, 1989, 224 с.

[31] Вощинин А. П. и др. Метод анализа данных при интервальной нестатистической ошибке. Завод. лаборатория 56, №7, 1990, 76-81.

[32] Белов В.М., Суханова В. А., Александрова С. Я. Определение констант скорости простых химических реакций методом центра неопределенности. Томск, 1988, Деп. в ОНИИТЭХИМ 17.06.88, №728-хп 88.

[33] ЭмАнуэль И. М., Кнорре Д. Г. Курс химической кинетики. Высш. шк., М., 1984.