CC BY
51
НАУЧНОЕ ИЗДАНИЕ МГТУ ИМ. Н. Э. БАУМАНА
НАУКА и ОБРАЗОВАНИЕ
Эл № ФС77 • 48211. Государственная регистрация №0421200025. ISSN 1994-0408
электронный научно-технический журнал
Новый конвективно-диффузионный метод глобальной
минимизации для решения обратных задач химической
кинетики
# 04, апрель 2013
Б01: 10.7463/0413.0569246
Федоров В. В.
УДК 519.6
Россия, Тольяттинский государственный университет
vvfmail@mail.ru
ВВЕДЕНИЕ
Обратная задача химической кинетики состоит в восстановлении по известной зависимости концентрации веществ от времени схемы реакции и констант скорости. Задача определения констант скоростей сводится к решению задачи минимизации отклонений между экспериментальными и расчетными данными [1].
Математические модели химической кинетики обычно представляют собой "жесткие" системы нелинейных дифференциальных уравнений, решение которых само по себе представляет сложную задачу. Как следствие, целевая функция может оказаться с "овражной" поверхностью и с несколькими локальными экстремумами.
Из-за больших вычислительных затрат, как в детерминированных так и в стохастических существующих методах, возникает необходимость применения алгоритмов с параллельными вычислениями [2-5] для чего, в свою очередь, требуется использование мощных многопроцессорных суперкомпьютеров [2], [3], [5].
Таким образом, вопрос о выборе эффективного метода глобальной минимизации многоэкстремальных функций нескольких переменных для идентификации математических моделей химической кинетики продолжает оставаться актуальным.
Предлагаемый метод основан на двух концепциях.
Первая концепция - сведение многомерной минимизации к одномерной минимизации, заложена в методе [6], где имитируется движение материальной точки по квазиоптимальным траекториям, рассчитываемых с помощью системы обыкновенных дифференциальных уравнений
с разными начальными условиями.
Модификацией последнего метода является так называемый метод тяжелого шарика, в котором моделируется движение шарика массой т по неровной поверхности с заданной начальной скоростью:
Рисунок 1. Движение тяжелого шарика по неровной поверхности
Здесь шарик с большой массой m по инерции проскакивает неглубокие минимумы. Однако при большой инерции он может проскочить и глобальный минимум, а при малой инерции может попасть в неглубокий локальный минимум. В результате требуется настройка нескольких параметров (m, v0, p) для конкретной целевой функции, что собственно представляет собой дополнительную задачу оптимизации.
Вторая концепция - это диффузионное сглаживание поверхности целевой функции, заложенная в методе DEM [7], [8], в котором целевая функция f(x) преобразуется в функцию F(x) без локальных экстремумов в результате решения уравнения диффузии
с начальными условиями:
Поиск глобального минимума ведется на сглаженной функции F (х, t) без
локальных минимумов, как изображено на рисунке 2. Но применение данного метода для решения обратных задач усложняется из-за вычисления вторых производных функционалов.
Рисунок 2. Исходная функция с локальными минимумами f(x*) и сглаженная функция
Обе эти концепции: сведение многомерной минимизации к одномерной и диффузионное сглаживание присутствуют в уравнении конвективной диффузии.
ОПИСАНИЕ МЕТОДА
Метод основан на интегрировании уравнений конвективной диффузии
с начальными условиями:
и граничными условиями:
(1)
(2)
(3)
до наступления стационарного режима.
В уравнениях переменные —искомые параметры; с10, Дс^ —величины для
определения области поиска минимума; коэффициенты диффузии, имеющие значения
0,01-ЮЛ; скорости конвективного движения, пропорциональные координатам
градиента целевой функции:
(4)где / - целевая функция.
Физическая интерпретация метода
В методе имитируется диффузионный перенос вещества с концентрацией с^ в многокомпонентном потоке, движущемся в гравитационном поле по неровной поверхности /(с1г... сп) (см. рисунок 3).
Рисунок 3. Физическая интерпретация метода
Распределение концентраций ¡-го вещества в движущемся потоке описывается уравнением конвективной диффузии:
Так как величины с1 с одной стороны являются концентрациями /-го вещества в
потоке, а с другой - координатами многомерного пространства, то векторы и дгас1с{ должны иметь вид:
(6)
(7)
Скорости же потоков зависят от направления и степени уклона поверхности, на основании чего их координаты будут:
или
(8)
Подставив выражения векторов в уравнение конвективной диффузии (5), получим уравнения вида (1).
Структура программы
Для реализации предлагаемого метода минимизации была разработана программа, состоящая из двух основных модулей - главного с процедурами расчета целевой функции и производных и модуля минимизации с процедурой расчета краевой задачи. Схематично структура программы представлена на рисунке .
Рисунок 4. Структура программы глобальной минимизации конвективно-диффузионным
методом
Влияние коэффициента диффузии
Дифференциальные уравнения необходимо интегрировать по возможности с малыми значениями коэффициентов диффузии Б. На рисунках 5, 6 представлены результаты минимизации функции
с разными коэффициентами Б.
Рисунок 5. Траектория конвективно-диффузионного потока по поверхности: а - при
Б=0,1; б - при Б=0,5
Рисунок 6. Распределение значений искомых переменных, функции и производных вдоль
траектории I: а - при Б=0,1; б - при Б=0,5
Оценка сложности метода минимизации
Сложность предлагаемого метода можно определить в виде:
где 1(п) — количество необходимых итераций, <р(п) —сложность вычисления
производной целевой функции, п — размерность минимизации, т — время одной итерации для одного уравнения конвективной диффузии.
При I = const и <р = const зависимость сложности от размерности (временная сложность) будет прямо пропорциональна n:
Для оценки временной сложности метода были выполнены тестовые расчеты по минимизации функции Розенброка
и сферической функции
с разными размерностями. Для сферической функции величины щ задавались
произвольно с помощью генератора случайных чисел (—3,8 < а{ < 3,8).
Интегрирование уравнений конвективной диффузии проводилось методом конечных разностей с количеством точек на отрезке равным 10 и коэффициенте диффузии 0,01 до наступления стационарного состояния с заданной точностью. Начальные и граничные условия были следующими:
< 1 Ci = 4 при I = -1; Q = -4 при I = 1.
Частные производные рассчитывались по формулам для функции Розенброка:
df (съ -сп+1) dci
400с,3 + 2Ci (1 - 200с!+1) - 2, i = 1 —200с1г_1 + 400с,3 + 2с, (101 - 200ci+1) - 2, 1 < i < п,
200ct - 200cf_lt i=n+ 1
для сферической функции:
-д-= - а,)
Результаты представлены в таблице 1 и в графическом виде на рисунке 7.
Таблица 1. Зависимость времени минимизации t от размерности п
Размерность, п Функция Розенброка Сферическая функция
Время, мс Целевая функция, f Время, мс Целевая функция, f
21 - - 172 4.6763 10-9
41 - - 218 1.201810-8
61 78 1.1606 10-1 296 2.5850 10-6
81 94 4.0052 10-2 343 2.3251 10-6
101 125 3.419510-2 421 4.5686 10-6
121 156 1.1805 10-2 468 1.7385 10-6
141 203 4.3424 10-2 561 2.1974 10-6
161 249 2.9390 10-2 749 2.8076 10-6
181 281 5.814310-2 718 4.5892 10-6
201 343 5.3230 10-3 889 2.8663 10-5
Рисунок 7. Зависимость времени минимизации t от размерности п t=f(n): слева - для функции Розенброка; справа - для сферической функции; числа в прямоугольниках - количество итераций.
Так как в данных примерах сложность вычисления производной целевой функции не зависит от п, а количество необходимых итераций практически не меняется, то
временная сложность действительно имеет порядок О (п).
К сожалению, в задачах параметрической идентификации математической модели химической кинетики сложность вычисления производной целевой функции
пропорциональна п:
так как количество дифференциальных уравнений и искомых констант скоростей химических реакций обычно бывают величинами одного порядка. Если подставить <р (п) в
уравнение (9), то временная сложность для таких задач будет порядка 0(п2). Тем не
менее, известно, что методы минимизации гладких выпуклых функций с применением производных, (например, градиентные методы) сходятся к минимуму со скоростью геометрической прогрессии.
Примеры тестирования многомерной минимизации предлагаемым методом приведены в [9]. Ниже приводится пример многопараметрической идентификации математической модели химической кинетики.
ПРИМЕР ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ МАТЕМАТИЧЕСКОЙ МОДЕЛИ ХИМИЧЕСКОЙ КИНЕТИКИ
В ряде работ, например в [10], [11], [12], для тестирования поведения жестких систем дифференциальных уравнений приводится система, взятая из уравнений химических реакций HIRES, предложенных Шефером [13] для объяснения "роста и дифференциации растительной ткани независимо от фотосинтеза при высоких уровнях светового облучения". Вопрос о численном решении жестких систем дифференциальных уравнений заслуживает отдельного внимания. Отметим лишь, что в данном случае система решалась по разностной схеме Розенброка с комплексными коэффициентами [14], [15], [16]. Интегрирование проводилось с шагом 0.1 на отрезке от 0 до 5 и с увеличенными шагами с количеством 100 на отрезке от 5 до 321.8122.
Математическая модель [12]:
с начальными условиями:
.
В работе [12] приведены значения концентраций С1И С4 при t = 1,2,3,4,5 и
значения всех концентраций в конце отрезка. Эти данные были приняты как экспериментальные. Остальные, принятые в качестве экспериментальных данных, были
рассчитаны с константами: к1 = 1.71; к2 = 0.43; к3 = 8.32; кА = 0.69; к5 = 0.035;
кв = 8.32; к7 = 280; к8 = 0.69; кэ = 0.69; к10 = 0.0007.
В обратной задаче требуется вычисление констант к^ по известным концентрациям. Целевая функция принята в виде:
где С^, — экспериментальные значения, С^ —рассчитанные значения /-го компонента на
7'-м шаге в процессе минимизации функционала, г^ — весовые коэффициенты.
Производные целевой функции вычислялись методом конечных разностей:
где 5 = 1" 10 4 (малое положительное число).
Для сведения многомерной минимизации целевой функции /(к1,...,к10) к
минимизации на одном отрезке решалась нестационарная краевая задача вида (1) с начальными условиями:
и граничными условиями:
где С; = 1пк,1,1 = 1,2, ...ДО.
Интегрирование уравнений проводилось методом конечных разностей с количеством точек на отрезке равным 16 и коэффициенте диффузии Б= 0,01 до значения
Г с 1 ■ Ю-4. На рисунке 8 изображены значения искомых параметров и целевой функции в одномерной области. Наименьшее значение f находится в точке 7.
/
Рисунок 8. Значения искомых параметров и целевой функции на отрезке
В таблице 2 представлены расчетные и фактические константы химических реакций.
Таблица 2. Расчетные и фактические значения параметров
Наименование Фактические значения Расчетные значения
к1 1.7100 1.7146
к2 0.4300 0.4277
кз 8.3200 8.7947
к4 0.6900 0.6923
к5 0.0350 0.0304
кб 8.3200 8.2662
к? 280.0000 293.0966
к8 0.6900 0.2428
к9 0.6900 0.9868
кю 0.0007 0.0023
На рисунках 9 , 10 изображены кинетические кривые, рассчитанные с фактическими и расчетными константами химических реакций.
Рисунок 9. Расчетные и "экспериментальные" кинетические кривые в начале отрезка
О
Ьы-1 ^
¿К-........ Д. ^_^^ ■ . ° П □ г, л ЛЧ
О 50 100 150 200 250 300
-С1 ■ експ
С2 о експ
С3х50 Д експ — С4 Т експ -СЕ + експ -06 □ експ -С 7x50 • експ —С8х100
Рисунок 10. Расчетные и "экспериментальные" кинетические кривые на всем отрезке
В таблице 3 редставлены расчетные, полученные в результате параметрической идентификации, и "экспериментальные", приведенные в [12].
Таблица 3. Расчетные и "экспериментальные" концентрации
Время, 1 Концентрации Значения
Расчетные "Экспериментальные"
1.0000 С1 2.558866 0.255493
С 4 0.458810 0.458520
2.0000 С1 0.123219 0.123520
С4 0.343855 0.344053
3.0000 С1 0.074510 0.074697
С 4 0.220323 0.220462
4.0000 С1 0.047866 0.047823
Q 0.139706 0.139692
5.0000 С1 0.031889 0.031652
С4 0.089957 0.089743
С1 0.000223 0.0007371
С2 0.000414 0.0001442
Сз 0.000158 0.0000589
321.8122 С4 0.000350 0.0011757
С5 0.002383 0.0023863
Сб 0.006276 0.0062390
Су 0.003103 0.0028500
С8 0.002597 0.0028500
ЗАКЛЮЧЕНИЕ
Разработан новый детерминированный метод глобальной минимизации многоэкстремальных функций нескольких переменных, основанный на решении нестационарного уравнения конвективной диффузии. Данным методом по разработанной программе решена задача параметрической идентификации математической модели химической кинетики. Результаты подтверждают применимость предлагаемого метода для решения подобных задач.
Список литературы
1. Полак Л.С., Гольденберг М.Я., Левицкий А.А. Вычислительные методы в химической кинетике. М.: Наука, 1984. 280 с.
2. Губайдуллин И.М., Рябов В.В., Тихонова М.В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислительные методы и программирование: новые вычислительные технологии. 2011. Т. 12, № 1. С. 137-145.
3. Линд Ю.Б. Применение суперкомпьютера для для решения обратных задач химической кинетики // Вычислительные технологии. 2006. Т. 11, Спец. вып. № 4 «Избранные доклады VIII Международного семинара-совещания "Кубатурные формулы и их приложения" (Улан-Удэ, август 2005 г.)». С. 76-80.
4. Завриев С.К., Перунова Ю.Н. Параллельные версии модифицированных методов покоординатного и градиентного спуска и их применение для решения некоторого класса задач глобальной оптимизации // Прикладная математика и информатика. М.: Диалог-МГУ. 2002. № 10.
5. Стронгин Р.Г., Гергель В.П., Баркалов К. А. Параллельные методы решения задач глобальной оптимизации // Изв. вузов. Приборостроение. 2009. Т. 52, № 10. С. 25-33.
6. Snyman A.J. A New and Dynamic Method for Unconstrained Optimization //Applied Mathematical Modelling. 1982. Vol. 6, no. 6. P. 449-462. http://dx.doi.org/10.1016/S0307-904X(82)80007-3
7. Piela L., Kostrowicki J., Scheraga H.A. The multiple-minima problem in the conformational analysis of molecules. Deformation of the potential energy hypersurface by the diffusion equation method // Journal of Physical Chemistry. 1989. Vol. 93. P. 3339-3346.
8. Hartman-Baker R.J. The Diffusion Equation Method for Global Optimization and its Application to Magnetotelluric Geoprospecting: Technical Report. 2005. Режим доступа: https://www.ideals.illinois.edu/handle/2142/11055 (дата обращения 09.01.2012).
9. Федоров В.В. Метод конвективно-диффузионной глобальной минимизации для многопараметрической идентификации математических моделей // Вектор науки Тольяттинского ГУ. 2012. № 3(21). С. 46-48.
10. Хайрер Э., Ваннер Г. Решение обыкновенных дифференциальных уравнений. Жесткие и дифференциально-алгебраические задачи: пер. с англ. М.: Мир. 1999. 685 с.
11. Холодов А.С., Лобанов А.И., Евдокимов А.В. Разностные схемы для решения жестких обыкновенных дифференциальных уравнений в пространстве неопределенных коэффициентов. Режим доступа:
http://crec.mipt.ru/study/materials/compmath/method/Zhestkie Syst.pdf (дата обращения 29.01.2013).
12. A.Emimal Kanaga Pushpam, D.Paul Dhayabaran. An Application of STWS Technique in Solving Stiff Non-linear System: 'High Irradiance Responses' (HIRES) of Photomorphogenesis // Recent Research in Science and Technology. 2011. Vol. 3 (6). P. 47-53
13. Sch'afer E. A new approach to explain the "High Irradiance Responses" of photomorphogenesis on the basis of phytochrome // J. of Math. Biology. 1975. Vol. 2. P. 41-56.
14. Днестровская Е.Ю., Калиткин Н.Н., Ритус И.В. Решение уравнений в частных производных схемами с комплексными коэффициентами // Математическое моделирование. 1991. Т. 3, № 9. С. 114-127.
15. Альшин А.Б., Альшина Е.А., Калиткин Н.Н., Корягина А.Б. Схемы Розенброка с комплексными коэффициентами для жестких и дифференциально-алгебраических систем // Журнал вычислительной математики и математической физики. 2006. Т. 46, № 8.
С. 1392-1414.
16. Rosenbrock H.H. Some general imlicit prosses for the numerical solution of differential equations // Computer jorn. 1963. Vol. 5. № 4. P. 329-330.
SCIENTIFIC PERIODICAL OF THE RAIJMAN MS TU
SCIENCE and EDUCATION
EL № FS77 - 48211. №0421200025. ISSN 1994-040S
electronic scientific and technical journal
A new convection-diffusion global minimization method for solving inverse problems of chemical kinetics # 04, April 2013 DOI: 10.7463/0413.0569246
Fedorov V.V.
Russia, Togliatti State University vvfmail@mail.ru
A new deterministic global minimization method based on the numerical solution of the convection-diffusion boundary problem was developed. With the help of this method, the inverse problem of chemical kinetics was solved. Problems of identification of speed constants in chemical reactions are usually based on minimization of functions. In most cases it is necessary to deal with "stiff systems of differential equations and, respectively, with many objective functions with ravine surface and local extremums, which makes application of existing multidimensional minimization methods complicated. Two main conceptions are included in this method: diffuse overcoming of local extremums and transforming multivariable minimization to one-dimensional minimization.
Publications with keywords: inverse problem, global minimization, parameters identification Publications with words: inverse problem, global minimization, parameters identification
References
1. Polak L.S., Gol'denberg M.Ia., Levitskii A.A. Vychislitel'nye metody v khimicheskoi kinetike [Computational methods in chemical kinetics]. Moscow, Nauka, 1984. 280 p.
2. Gubaidullin I.M., Riabov V.V., Tikhonova M.V. Primenenie indeksnogo metoda global'noi optimizatsii pri reshenii obratnykh zadach khimicheskoi kinetiki [Application of the global optimization index method to solving inverse problems of chemical kinetics]. Vychislitel'nye metody i programmirovanie: novye vychislitel'nye tekhnologii [Numerical methods and programming. Advanced Computing], 2011, vol. 12, no. 1, pp. 137-145.
3. Lind Iu.B. Primenenie superkomp'iutera dlia dlia resheniia obratnykh zadach khimicheskoi kinetiki [Application of supercomputer to solution of inverse problems of chemical kinetics]. Vychislitel'nye tekhnologii [Computational Technologies], 2006, vol. 11, spec. iss. no. 4, pp. 76-80.
4. Zavriev S.K., Perunova Iu.N. Parallel'nye versii modifitsirovannykh metodov pokoordinatnogo i gradientnogo spuska i ikh primenenie dlia resheniia nekotorogo klassa zadach
global'noi optimizatsii [Parallel versions of the modified coordinate and gradient descent methods and their application to a class of global optimization problems]. Prikladnaia matematika i informatika, 2002, no. 10. (Trans. version: Computational Mathematics and Modeling, 2003, vol. 14, no. 2, pp. 108-122. DOI: 10.1023/A:1022951305895).
5. Strongin R.G., Gergel' V.P., Barkalov K. A. Parallel'nye metody resheniia zadach global'noi optimizatsii [Parallel methods for global optimization problem solvin]. Izv. vuzov. Priborostroenie, 2009, vol. 52, no. 10, pp. 25-33.
6. Snyman A.J. A New and Dynamic Method for Unconstrained Optimization. Applied Mathematical Modelling, 1982, vol. 6, no. 6, pp. 449-462. http://dx.doi.org/10.1016/S0307-904X(82)80007-3
7. Piela L., Kostrowicki J., Scheraga H.A. The multiple-minima problem in the conformational analysis of molecules. Deformation of the potential energy hypersurface by the diffusion equation method. Journal of Physical Chemistry, 1989, vol. 93, pp. 3339-3346.
8. Hartman-Baker R.J. The Diffusion Equation Method for Global Optimization and its Application to Magnetotelluric Geoprospecting: Technical Report. 2005. Available at: https://www.ideals.illinois.edu/handle/2142/11055 , accessed 09.01.2012.
9. Fedorov V.V. Metod konvektivno-diffuzionnoi global'noi minimizatsii dlia mnogoparametricheskoi identifikatsii matematicheskikh modelei [A new convection-diffuzion global method of minimization for multivariable identification of the mathematical models]. Vektor nauki Tol'iatti SU, 2012, no. 3 (21), pp. 46-48.
10. Hairer E., Wanner G. Solving Ordinary Differential Equations 2: Stiff and Differential-Algebraic Problems. Berlin, Springer-Verlag, 1991. 601 p. (Springer Series in Computational Mathematic). (Russ. ed.: Khairer E., Vanner G. Reshenie obyknovennykh differentsial'nykh uravnenii. Zhestkie i differentsial'no-algebraicheskie zadachi. Moscow, Mir, 1999. 685 p.).
11. Kholodov A.S., Lobanov A.I., Evdokimov A.V. Raznostnye skhemy dlia resheniia zhestkikh obyknovennykh differentsial'nykh uravnenii v prostranstve neopredelennykh koeffitsientov. Available at:
http://crec.mipt.ru/study/materials/compmath/method/Zhestkie_Syst.pdf , accessed 29.01.2013.
12. A.Emimal Kanaga Pushpam, D.Paul Dhayabaran. An Application of STWS Technique in Solving Stiff Non-linear System: 'High Irradiance Responses' (HIRES) of Photomorphogenesis. Recent Research in Science and Technology, 2011, vol. 3 (6), pp. 47-53
13. Sch'afer E. A new approach to explain the "High Irradiance Responses" of photomor-phogenesis on the basis of phytochrome. J. of Math. Biology, 1975, vol. 2, pp. 41-56.
14. Dnestrovskaia E.Iu., Kalitkin N.N., Ritus I.V. Reshenie uravnenii v chastnykh proizvodnykh skhemami s kompleksnymi koeffitsientami [Solving of partial differential equations by schemes with complex coefficients]. Matematicheskoe modelirovanie, 1991, vol. 3, no. 9, pp. 114-127.
15. Al'shin A.B., Al'shina E.A., Kalitkin N.N., Koriagina A.B. Skhemy Rozenbroka s kompleksnymi koeffitsientami dlia zhestkikh i differentsial'no-algebraicheskikh sistem [Rosenbrock schemes with complex coefficients for stiff and differential algebraic systems]. Zhurnal vychislitel'noi matematikii matematicheskoi fiziki, 2006, vol. 46, no. 8, pp. 1392-1414. (Trans. version: Computational Mathematics and Mathematical Physics, 2006, vol. 46, no. 8, pp. 1320-1340. DOI: 10.1134/S0965542506080057 ).
16. Rosenbrock H.H. Some general imlicit prosses for the numerical solution of differential equations. Computer jorn., 1963, vol. 5, no. 4, pp. 329-330.
