Другие журналы

научное издание МГТУ им. Н.Э. Баумана

НАУКА и ОБРАЗОВАНИЕ

Издатель ФГБОУ ВПО "МГТУ им. Н.Э. Баумана". Эл № ФС 77 - 48211.  ISSN 1994-0408

Новый конвективно-диффузионный метод глобальной минимизации для решения обратных задач химической кинетики

# 04, апрель 2013
DOI: 10.7463/0413.0569246
Файл статьи: Федоров_P.pdf (791.82Кб)
автор: Федоров В. В.

УДК 519.6

Россия, Тольяттинский государственный университет

vvfmail@mail.ru

 

ВВЕДЕНИЕ

 

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

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

Из-за больших вычислительных затрат, как в детерминированных так и в стохастических существующих методах, возникает необходимость применения алгоритмов с параллельными вычислениями [2-5] для чего,  в свою очередь, требуется использование мощных многопроцессорных суперкомпьютеров [2], [3], [5].

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

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

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

с разными начальными условиями.

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

 

Рисунок 1. Движение тяжелого шарика по неровной поверхности

 

Здесь шарик с большой массой m по инерции проскакивает неглубокие минимумы. Однако при большой инерции он может проскочить и глобальный минимум, а при малой инерции может попасть в неглубокий локальный минимум. В результате требуется настройка нескольких параметров (m, v0, p) для конкретной целевой функции, что собственно представляет собой дополнительную задачу оптимизации.

Вторая концепция – это диффузионное сглаживание поверхности целевой функции, заложенная в методе DEM [7], [8], в котором целевая функция f(x) преобразуется в функцию F(x) без локальных экстремумов в результате решения уравнения диффузии

с начальными условиями:

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

 

 

Рисунок 2. Исходная функция с локальными минимумами  и сглаженная функция

Обе эти концепции: сведение многомерной минимизации к одномерной и диффузионное сглаживание присутствуют  в уравнении конвективной диффузии.

 

ОПИСАНИЕ МЕТОДА

Метод основан на интегрировании уравнений конвективной диффузии

                       (1)

с начальными условиями:

                       (2)

и граничными условиями:

       (3)

до наступления стационарного режима.

В уравнениях  переменные искомые параметры;  величины для определения области поиска минимума; – коэффициенты диффузии, имеющие значения 0,01÷0,1; – скорости конвективного движения, пропорциональные координатам градиента целевой функции:

                                  (4)

где f – целевая функция.

Физическая интерпретация метода

В методе имитируется диффузионный перенос вещества с концентрацией  в многокомпонентном потоке, движущемся в гравитационном поле по неровной поверхности  (см. рисунок 3).

Рисунок2-3.png

 

Рисунок 3. Физическая интерпретация метода

 

Распределение концентраций i-го вещества в движущемся потоке описывается уравнением конвективной диффузии:

               (5)

Так как величины  с одной стороны являются концентрациями i-го вещества в потоке, а с другой – координатами многомерного пространства, то векторы  и  должны иметь вид:

                                             (6)

                                 (7)

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

или

      (8)

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

Структура программы

Для реализации предлагаемого метода минимизации была разработана программа, состоящая из двух основных модулей – главного с процедурами расчета целевой функции и производных и модуля минимизации с процедурой расчета краевой задачи. Схематично структура программы представлена на рисунке 4.

 

 

Рисунок 4. Структура программы глобальной минимизации конвективно-диффузионным методом

 

Влияние коэффициента диффузии

Дифференциальные уравнения  необходимо интегрировать по возможности  с  малыми значениями коэффициентов диффузииD. На рисунках 5, 6 представлены результаты минимизации функции

с разными коэффициентами D.

 

 

Рисунок 5. Траектория конвективно-диффузионного потока по поверхности: а – при D=0,1; б – при D=0,5

 

 

Рисунок 6. Распределение значений искомых переменных, функции и  производных вдоль траектории l: а – при D=0,1; б – при D=0,5

 

Оценка сложности метода минимизации

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

         (9)

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

При  и  зависимость сложности от размерности (временная сложность) будет прямо пропорциональна n:

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

 и сферической функции

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

Интегрирование уравнений конвективной диффузии проводилось методом конечных разностей с количеством точек на отрезке равным 10 и коэффициенте диффузии D= 0,01 до наступления стационарного состояния с заданной точностью.

Начальные и граничные условия были следующими:

Частные производные рассчитывались по формулам

для функции Розенброка:

для сферической функции:

 

Результаты представлены в таблице 1 и в графическом виде на рисунке 7.

Таблица 1. Зависимость времени минимизации tот размерности n

Размерность, n

ФункцияРозенброка

Сферическая функция

Время, мс

Целевая функция, f

Время, мс

Целевая функция, f

21

-

-

172

4.6763.10-9

41

-

-

218

1.2018.10-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.4195.10-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.8143.10-2

718

4.5892.10-6

201

343

5.3230.10-3

889

2.8663.10-5

 

time.emftime.emf

 

Рисунок 7. Зависимость времени минимизации tот размерности n t=f(n):

слева – для функции Розенброка; справа – для сферической функции; числа в прямоугольниках – количество итераций.

 

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

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

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

            Примеры тестирования многомерной минимизации предлагаемым методом приведены в [9]. Ниже приводится пример многопараметрической идентификации математической модели химической кинетики.

 

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

В ряде работ, например в [10], [11], [12], для тестирования поведения жестких систем дифференциальных уравнений приводится система, взятая из уравнений химических реакций HIRES, предложенных Шефером [13] для объяснения “роста и дифференциации растительной ткани независимо от фотосинтеза при высоких уровнях светового облучения”. Вопрос о численном решении жестких систем дифференциальных уравнений заслуживает отдельного внимания. Отметим лишь, что в данном случае  система решалась по разностной схеме Розенброка с комплексными коэффициентами [14], [15], [16]. Интегрирование проводилось с шагом 0.1 на отрезке от 0 до 5 и с увеличенными шагами с количеством 100 на отрезке от 5 до 321.8122.

Математическая модель [12]:

с начальными условиями: .

В работе [12] приведены значения концентраций и  при  и значения всех концентраций в конце отрезка. Эти данные были приняты как экспериментальные. Остальные, принятые в качестве экспериментальных данных, были рассчитаны с константами: ; ; ; ; ; ; ; ; ; .

В обратной задаче требуется вычисление констант  по известным концентрациям.

Целевая функция принята в виде:

где  экспериментальные значения,   рассчитанные значения i-го компонента на j-м шаге в процессе  минимизации функционала, .

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

где  (малое положительное число).

Для сведения многомерной минимизации целевой функции  к минимизации на одном отрезке решалась нестационарная краевая задача вида (1) с начальными условиями:

и граничными условиями:

где

Интегрирование уравнений  проводилось методом конечных разностей с количеством точек на отрезке равным 16 и коэффициенте диффузии D= 0,01 до значения . На рисунке 8  изображены значения искомых параметров и целевой функции в одномерной области. Наименьшее значение  находится в точке 7.

 

 

Рисунок 8. Значения искомых параметров и целевой функции на отрезке

 

В таблице 2 представлены расчетные и фактические константы химических реакций.

            Таблица 2. Расчетные и фактические значения параметров

Наименование

Фактические значения

Расчетные значения

k1

1.7100

  1.7146

k2

0.4300

  0.4277

k3

8.3200

  8.7947

k4

0.6900

  0.6923

k5

0.0350

  0.0304

k6

8.3200

  8.2662

k7

280.0000

293.0966

k8

0.6900

  0.2428

k9

0.6900

  0.9868

k10

0.0007

  0.0023

 

На рисунках 9 , 10 изображены кинетические кривые, рассчитанные с фактическими и расчетными константами химических реакций.

 

Рисунок 9. Расчетные и “экспериментальные” кинетические кривые в начале отрезка

 

 

Рисунок 10. Расчетные и “экспериментальные” кинетические кривые на всем отрезке

 

В таблице 3 редставлены расчетные, полученные в результате параметрической идентификации, и “экспериментальные”, приведенные в [12].

              Таблица 3. Расчетные и “экспериментальные” концентрации

Время, t

Концентрации

Значения

Расчетные

“Экспериментальные”

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

С4

0.139706

0.139692

5.0000

С1

0.031889

0.031652

С4

0.089957

0.089743

321.8122

С1

0.000223

0.0007371

С2

0.000414

0.0001442

С3

0.000158

0.0000589

С4

0.000350

0.0011757

С5

0.002383

0.0023863

С6

0.006276

0.0062390

С7

0.003103

0.0028500

С8

0.002597

0.0028500

 

ЗАКЛЮЧЕНИЕ

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

 

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

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

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 photomor-phogenesis 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.

Поделиться:
 
ПОИСК
 
elibrary crossref ulrichsweb neicon rusycon
 
ЮБИЛЕИ
ФОТОРЕПОРТАЖИ
 
СОБЫТИЯ
 
НОВОСТНАЯ ЛЕНТА



Авторы
Пресс-релизы
Библиотека
Конференции
Выставки
О проекте
Rambler's Top100
Телефон: +7 (915) 336-07-65 (строго: среда; пятница c 11-00 до 17-00)
  RSS
© 2003-2024 «Наука и образование»
Перепечатка материалов журнала без согласования с редакцией запрещена
 Тел.: +7 (915) 336-07-65 (строго: среда; пятница c 11-00 до 17-00)