CC BY
155
УДК 021.8 + 025.1 ББК 78.34
РЕШЕНИЕ ОПТИМИЗАЦИОННЫХ ЗАДАЧ ПРИ ХИМИЧЕСКИХ ПРЕВРАЩЕНИЯХ
Майков И. Л.1, Зайченко В. М.2
(ФГБУН Объединенный Институт высоких температур РАН, Москва)
Разработан алгоритм решения оптимизационных задач определения параметров химических превращений на основе метода Хука и Дживса, использующий одномерную минимизацию вдоль координатных направлений. Алгоритм включает два типа поиска - исследующий поиск и поиск по образцу с введением ограничений. Одномерная оптимизация проводится на основе алгоритма Брента. Для вычисления целевой функции решается система жестких дифференциальных уравнений с использованием алгоритм DIFSUB.
Ключевые слова: математическое моделирование, оптимизация, химическая кинетика, нелинейное программирование, методы поиска.
1. Введение
Необходимость вовлечения в энергобаланс местных топливно-энергетических ресурсов (торфа, низкосортных углей, растительных сельскохозяйственных отходов, отходов лесозаготовительной и деревообрабатывающей промышленности) сти-
1 Игорь Леонидович Майков, доктор физико-математических наук, в.н.с. ([email protected]).
2 Виктор Михайлович Зайченко, доктор технических наук, зав.лаб. ^аН^@оМгап. ги).
мулирует разработку установок и схем энергетического использования биоресурсов. Одно из наиболее перспективных направлений - термохимическая переработка растительного сырья в кондиционный энергетический газ (процессы газификации и пиролиза) [9]. Структурная схема такого энерготехнологического комплекса представлена на рис. 1. Добавление в схему блока торификации ^ог^ас^оп) позволяет наряду с основными продуктами комплекса (электроэнергия и тепло) получать Ьу-продукт - топливные пеллеты с улучшенными физикомеханическими свойствами [11].
Сырье ■ >
1
1
з
к.
4
^
Электричество
Тепло
=>
Товарный продукт =>
Рис. 1. Структурная схема энерготехнологического комплекса на растительных отходах: 1 — термохимический реактор;
2 — газопоршневая мини-ТЭЦ; 3 — котел с кипящим слоем;
4 — блок торификации
Очевидно, что оптимальный состав и характеристики компонентов комплекса зависят от многих факторов (в первую очередь от сырьевой базы) и могут быть достигнуты только в результате решения задачи многопараметрической оптимизации (примеры оптимизации и системы управления энергетического комплекса на основе газопоршневой мини-ТЭЦ приведены в [4]). Наряду с общей задачей оптимизации энерготехнологического комплекса возникают задачи оптимизации и управления отдельными блоками, основным из которых является термохимический реактор. Характерной особенностью термохимического реактора является наличие возмущающих параметров случайного характера, а именно, неизвестен точный компонентный состав сырья и имеются только некие усредненные характери-
стики, которые варьируют в широком диапазоне даже для определенного вида сырья. Кроме того, имеется тенденция изменения свойств во времени (например, торф при хранении изменяет свои свойства). В настоящее время существует большое количество работ, посвященных математическому моделированию термического разложения биомассы, в которых рассматриваются схемы превращения с различной степенью детализации (количества химических реакций). В результате решения обратной задачи (задачи оптимизации) по экспериментальным данным получают основные параметры процесса разложения конкретного вида биомассы. К сожалению, эти параметры не являются универсальными, и использование их для других видов биомассы и других режимов (например, темп нагрева) может привести к серьезным ошибкам в определении, например, времени процесса. Соответственно, построенная на таких параметрах система управления будет по крайней мере неэффективной.
Альтернативный способ построения системы управления состоит в том, что параметры химических реакций не являются постоянными величинами, а изменяются во времени (по длине реактора) согласно случайному характеру изменения состава сырья и, следовательно, протеканию процесса. То есть необходимо в каждый момент времени иметь распределение по длине реактора, например, оставшейся (непрореагировавшей) биомассы. Тогда, решая обратную задачу (оптимизационную задачу), можно получить параметры химических превращений, описывающих процесс в этот момент времени. И далее определить оптимальное воздействие на эту систему с целью перевода процесса в оптимальный режим (по выбранному критерию). Таким образом, алгоритм решения оптимизационной задачи для определения параметров процесса термического разложения нужно рассматривать как часть алгоритмов, включаемых в систему управления термохимического реактора, что накладывает жесткие условия на его быстродействие в связи с необходимостью многократного и быстрого решения локальных задач оптимизации.
Целью настоящей работы является построение эффективного алгоритма решения оптимизационной задачи для определения параметров термической переработки. В работе рассмотрены алгоритмы поиска [1, 9] (метод циклического покоординатного спуска и его модификации, метод Хука и Дживса, метод Розенброка). Выбор этих алгоритмов в качестве базовых для исследования продиктован возможностью включения в алгоритмы ограничений, а также особенностями построения этих алгоритмов, характерными для большинства алгоритмов поиска. Проведены численные расчеты типичной задачи оптимизации процесса термической деструкции.
2. Постановка задачи
При количественном изучении процесса термической деструкции приходится решать следующую задачу оптимизации [5, 6]: известна схема превращений и экспериментальные кривые компонентов; при этом константы скоростей всех стадий неизвестны. Требуется найти численные значения констант таким образом, чтобы они приводили к удовлетворительному описанию эксперимента.
Подобные задачи оптимизации, сформулированные математически, могут быть объединены под общим названием задача нелинейного программирования [8]. Методы решения таких задач весьма разнообразны и не все из них нашли применение для исследования процессов термического разложения в связи со спецификой задачи.
Рассмотрим кратко математическую модель разложения вещества по N каналам. В предположении экспоненциальной зависимости скорости реакций от температуры [5] разложение по i-му каналу можно представить в виде [2, 3]
dzt , ( О
—- = -ki exp-----—
(1) dt г Я T (t)
Z, (0) = 1.
Шг
где / = 1, ..., N Zi - концентрация /-го компонента; ^ - время; ДО - заданная функция температуры от времени; ^-, Е/, т/ -параметры оптимизации (для /-го канала разложения ^ - предэкспоненциальный множитель; Е/ - энергия активации; т/ - порядок реакции). Дополнительно вводится параметр оптимизации п/ - доля /-го компонента. Из физических соображений на параметры оптимизации накладываются следующие ограничения:
(2) 0 < п/ < 1, Ь > 0, Е/ > 0, т/ > 1.
Экспериментальная кинетическая кривая имеет вид
(3) у(^) = Уу,
где у - у-я экспериментальная точка.
Тогда функцию ошибок можно представить в виде
М ( N
(4) / = ЁI Уу , Е/, т)
У=1V /=1
где М - количество экспериментальных точек.
Введем вектор х(п/, Ь, Е/, тг) = х(хь х2, ..., х„). Сформулируем следующую оптимизационную задачу:
найти такие х, при которых целевая функция (4) при ограничениях (2) достигает минимума
(5) / (х) ^ тт.
3. Предварительный анализ
В общем случае функция / зависит от 4N - 1 переменных
N
(так как ^, где ц0 - доля разлагающихся компонент,
/=1
является характеристикой биомассы и определяется экспериментально). Как видно из постановки задачи (5), для вычисления функции / при заданном векторе х(п/, Ь, Е/, тг) требуется решение задачи Коши (1) для системы N уравнений ^ определяется из решения дифференциального уравнения (1)). Система (1) представляет собой жесткую систему уравнений [7], что физически соответствует разложению по разным каналам с существенно различными Ег (энергия активации). То есть для решения
задачи (5) необходимо дополнительно иметь эффективный алгоритм решения системы жестких уравнений. В работе использовался линейный многошаговый метод с автоматическим выбором шага, реализованным в алгоритме DIFSUB [10].
По существу методы поиска заключаются в следующем: при заданном векторе х определяется допустимое направление d и функция / минимизируется вдоль направления d одним из методов одномерной оптимизации. Применение методов одномерной оптимизации в методе поиска позволяет ввести ограничения на переменные для решения задачи (5) условной оптимизации. В работе использовался алгоритм Брента [7], в основе которого лежит комбинация методов золотого сечения и последовательной параболической интерполяции.
4. Алгоритмы методов поиска
Метод циклического покоординатного спуска. В этом методе в качестве направлений поиска используются координатные векторы. Метод осуществляет поиск вдоль направлений dl, ..., d„, где dj - вектор, все компоненты которого, за исключением у-й, равны нулю. При поиске по направлению dj меняется только переменная х, в то время как все остальные переменные остаются зафиксированными.
Общая схема метода:
1. Определение в качестве d1, ..., d„ координатных направлений. Выбор нулевого приближения х® (к = 0). Задание у1 = х®.
2. Решение оптимизационной задачи /(у, + в dj) для определения ву. Присвоение у;+1 = у;- + в dj. Если у < п, то у = у + 1 и переход к п. 2, если у = п, то переход к п. 3.
3. Задание хда = уп+1. Проверка условия ]Дх(М)) -/хда)| < е. Если условие выполняется - выход, в противном случае -k = k + 1 и переход к п. 2
Метод циклического покоординатного спуска с ускоряющим шагом. Метод циклического покоординатного спуска при минимизации дифференцируемой функции сходится к точке с нулевым значением градиента. В отсутствие дифференцируемо-
сти метод может остановиться в неоптимальной точке [8], что объясняется наличием оврага, вызванного недифференцируемостью. Эта трудность может быть преодолена поиском вдоль
1 Ck+1) (с) ^+1) ^+1) I р Л
направления d = х1 ' - х1' и вычисления х1 ' = х1 ' + в d, где
в - ускоряющий коэффициент. Обычно эмпирическим путем устанавливается, что такой шаг делается на каждой р-й итерации. Такая модификация метода циклического покоординатного спуска может ускорять сходимость, в частности, когда последовательность точек образует зигзагообразную траекторию вдоль дна оврага.
Метод Хука и Дживса. Метод Хука и Дживса [8] осуществляет два типа поиска - исследующий поиск и поиск по образцу.
Общая схема метода:
1. Определение в качестве d1, ..., координатных направлений. Выбор нулевого приближения х® (к = 0). Задание у1 = х®.
2. Исследующий поиск. Решение оптимизационной задачи /(у, + в dj) для определения ву. Присвоение у;+1 = у, + ву dj. Если у < п, то у = у + 1 и переход к п. 2, если у = п, то переход к п. 3.
3. Поиск по образцу. Задание х1^1 = у„+1. Проверка условия JДxCk+1)) - .Дхда)| < е. Если условие выполняется - выход, если нет, то задать d = х1^1 - х®. Решение оптимизационной задачи /У'М'1 + АД) для определения X и вычисления у1 = х1^1 + ?id, присвоение k = k + 1 и переход к п. 2.
Метод Розенброка. На каждой итерации процедура осуществляет итеративный поиск вдоль п линейно независимых и ортогональных направлений. Когда получена новая точка в конце итерации, строится новое множество ортогональных векторов с использованием процедуры Грамма-Шмидта [8].
1. Определение в качестве Д1, ..., координатных направлений. Выбор нулевого приближения х® (к = 0). Задание у1 = х®.
2. Решение оптимизационной задачи /(у, + в Д,) для определения ву. Присвоение у;+1 = у, + ву Ду. Если , < п, то , =, + 1 и переход к п. 2, если, = п, то переход к п. 3.
3. Задание х|'М) = уп+1. Проверка условия |Дх|'М)) -/хда)| < е. Если условие выполняется - выход, если нет, то построение
нового множества линеино независимых и взаимно ортогональ-
^ 1 пв'М
ных направлении а;- в соответствии с правилами d ,, если р = 0,
(6)
£ рД., если рJ * 0,
. *=/
а, если j = 1,
j 7 ^ У
а,-£ (а}аг )аг, если / > 2,
ь,.
и переход к п. 2 с использованием в качестве а, нового базиса
а/т
Замечание. Во всех методах (кроме метода циклического покоординатного спуска) существует возможность выхода параметров за наложенные ограничения при минимизации вдоль какого-то направления (например, поиск по образцу в методе Хука и Дживса или последовательная минимизация вдоль /-го направления в методе Розенброка).
Рассмотрим в качестве примера введение ограничении в методе Хука и Дживса. При решении оптимизационной задачи дх(£+1) + ^а) при поиске по образцу для определения X необходимо определить ограничения на X (Хт/п и Хтах) через заданные ограничения на параметры хт/п(а1, а2, ..,, ап). Это можно сделать, используя следующую процедуру
1. Определение вектора А = х(4г+1) - хда = (Дь Д2, ..., Дп).
2. Присвоение / = 1, Хт/п = 1010, Хтах = -1010
3. Если Д, < 0, то я =
Ш1П
^ х+1 -аі abs(Дi) ’
Я
если Д, < 0, то
= шах
- хк+1
Я
abs(Д/.) ’
4. Если / = п, то выход, в противном случае - / = / + 1 и переход к п. 3.
Ь
,=1
ь
а
5. Результаты и обсуждение
В качестве исходных данных для решения задачи оптимизации использовались экспериментальные данные по разложению измельченной древесины березы [3], (количество экспериментальных точек М = 177). Предполагалась четырехканальная схема разложения (т.е. N = 4, количество параметров оптимизации п = 15, величина По равна 0,822). В качестве начального приближения использовалось решение оптимизационной задачи по одному каналу, цг = 0,205, i = 1, ..., 4.
Для различных методов оптимизации определялось время достижения критерия выхода
М ( 4 ^2
(7) I = 1 УJ-Ът?, {к, Е,, т,)
У=1
V 1=1
< 8 =10-3.
Результаты расчетов с использованием различных методов поиска приведены на рис. 2, время расчетов нормировалось на время расчета по методу циклического покоординатного спуска.
Рис. 2. Безразмерное время расчета оптимизационной задачи: 1 - метод циклического покоординатного спуска; 2-5 - метод циклического покоординатного спуска с ускоряющим шагом (2 - в = 0,5; 3 - в = 1; 4 - в = 2,5; 5 - в = 5);
6 - метод Хука и Дживса
Время расчета по методу циклического покоординатного спуска составляло 22 мин. Расчеты проводились на компьютере с двухядерным процессором Core 2 Duo с тактовой частотой 2,4 ГГц и производительностью 19,2 гигафлопс.
Результаты расчета параметров методом Хука и Дживса представлены в таблице 1.
Таблица 1. Параметры оптимизации
i ki Ei mi Пі
1 1,525 • 1011 18250,8 2,5538 0,2919
2 3,968 • 1026 41540,7 1,0002 0,2773
3 3,857 • 1016 23831,1 1,5275 0,1305
4 0,0399 2543,9 1,0002 0,1223
В таблице 2 приведены максимальные абсолютные относительные ошибки параметров оптимизации (%), полученных различными алгоритмами по сравнению с алгоритмом Хука и Дживса. Максимальное расхождение не превышает 1,5%.
Таблица 2. Максимальные абсолютные относительные ошибки
%
i ki Ei mi Пі
1 1,48 0,35 1,14 0,23
2 1,43 0,22 3,99 • 10-3 0,46
3 1,11 0,43 1,42 0,73
4 1,47 0,57 3,99 • 10-3 0,32
Метод циклического покоординатного спуска с ускоряющим шагом уменьшает время счета. В расчетах ускоряющий шаг использовался после каждой итерации. Но выбор параметра ускорения остается произвольным, т.е. для каждой конкретной экспериментальной кривой будет свой параметр ускорения, устанавливаемый эмпирически. Метод Розенброка не дал удовлетворительных результатов. Основная причина этого состоит в том, что при построении системы ортогональных векторов с использованием процедуры Грамма-Шмидта (6) возникают и
накапливаются ошибки из-за существенной разницы величин параметров оптимизации (см. таблицу 1 -параметры различаются на 28 порядков). По-видимому, аналогичное поведение будет наблюдаться и для других алгоритмов прямого поиска (например, метод деформируемого многогранника [8]), в которых необходимо выполнение операций над векторами (вычисление центра тяжести в методе деформируемого многогранника).
Наименьшее время достижения оптимальных параметров дает метод Хука и Дживса. Это связано с тем, что в исследующем поиске выполняется циклический поиск вдоль направлений координатных векторов (величины параметров оптимизации не оказывают влияния на сходимость метода), а поиск по образцу снимает взаимозависимость переменных.
В работе дополнительно исследовались вышеперечисленные прямые методы с использованием постоянных шагов по направлениям, вместо одномерной минимизации. Такие методы оказались неэффективными (временя расчета на порядки больше). Это опять же связано с огромной разницей в величинах параметров оптимизации, так как априори неизвестно, какие из параметров будут иметь наибольшие значения, какие - наименьшие, и вопрос выбора величины шага для конкретного параметра остается открытым.
6. Заключение
В работе показано, что разработанный алгоритм решения оптимизационных задач термической деструкции на основе метода Хука и Дживса обладает наибольшей эффективностью по сравнению с другими методами поиска. Это связано в первую очередь с тем, что параметры оптимизации имеют величины, отличающиеся на несколько десятков порядков, что приводит к значительным ошибкам при реализации других алгоритмов поиска.
Представленный алгоритм позволяет многократного и быстро решать локальные задачи оптимизации и может быть
использован как часть алгоритма системы управления термохимическим реактором.
Работа выполнялась в рамках контракта с Министерством
образования и науки РФ № 16.526.11.6017.
Литература
1. БАЗАРА М., ШЕТТИ К. Нелинейное программирование. Теория и алгоритмы. Пер. с англ. - М.: Мир, 1982. - 583 с.
2. ДИРЕКТОР Л.Б., ЗАЙЧЕНКО В.М., КОСОВ В.Ф., МАЙКОВ И.Л., СИНЕЛЬЩИКОВ В.А., СОКОЛ Г.Ф. Определение константы скорости образования пироуглерода в потоке бутана // Известия Академии наук. Энергетика. -2009. - №3. - С. 79-87.
3. МАЙКОВ И.Л., СИНЕЛЬЩИКОВ В.А., ФЕДЮХИН В.Ф. Исследование термического распада органического сырья растительного происхождения // Труды пятой российской национальной конференции по теплообмену, Москва, 25-29 октября 2010 г. Том 3. - С. 262-264.
4. МАЙКОВ И.Л., ДИРЕКТОР Л.Б. Решение задач оптимизации и управления гибридными энергетическими комплексами в структуре распределенной генерации // Управление большими системами. - 2011. - №35. - С. 250-264.
5. ПОЛАК Л.С. ГОЛЬДЕНБЕРГ М.Я. ЛЕВИЦКИЙ А.А. Вычислительные методы в химической кинетике. - М.: Наука, 1979. - 280 с.
6. Применение вычислительной математики в химической и физической кинетике / Под редакцией Л.С. Полака. - М.: Наука, 1969. - 279 с.
7. ФОРСАЙТ Д., МАЛЬКОЛЬМ М., МОУЛЕР К. Машинные методы математических вычислений. - М.: Мир, 1980. -280 с.
8. ХИММЕЛЬБЛАУ Д. Прикладное нелинейное программирование. Пер. с англ. - М.: Мир, 1975. - 536 с.
9. BASU P., KAUSHAL P. Modeling of Pyrolysis and Gasification of Biomass in Fluidized Beds: A Review // The Berkeley Elec-
tronic Press. Chemical Product and Process Modeling. - 2009. -Vol. 4. Issue 1., Article 21. - P. 137-172.
10. GEER C.W. DIFSUB for solution of ordinary differential equation // Communication of the ACM. - 1971. - Vol. 14, №3. -P. 341-354.
11. ZAICHENKO V.M., KOSOV V.V., KOSOV V.F., SINELSCHIKOV VA. Torrefaction and synthesis gas production // The Proceedings of 19th European Biomass Conference and Exhibition. 6-11 June 2011, Berlin, Germany. -P.2011-2014.
SOLVING OPTIMIZATION PROBLEMS IN CHEMICAL REACTIONS
Igor Maikov, Joint Institute for High Temperatures of RAS, Moscow, professor ([email protected]).
Victor Zaitchenko, Joint Institute for High Temperatures of RAS, Moscow, professor ([email protected]).
Abstract: An algorithm for solving optimization problems of chemical kinetics on the basis of Hooke and Jeeves, which uses a one-dimensional minimization along coordinate directions. The algorithm consists of two types of search - exploring search and pattern-matching search with restrictions. The one-dimensional optimization algorithm is based on the Brent’s algorithm. To calculate the objective function a system of stiff differential equations is solved with DIFSUB algorithm.
Keywords: mathematical modeling, optimization, chemical
kinetics, nonlinear programming, methods of retrieval.
Статья представлена к публикации членом редакционной коллегии Н. Н. Бахтадзе
