Двухкритериальная идентификация кинетических параметров реакции гидроалюминирования олефинов алкилаланами Текст научной статьи по специальности «Математика»

научное издание мгту им. н. э. баумана

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

Эл № ФС77 • 48211. Государственная регистрация №0421200025. ISSN 1994-0408

электронный научно-технический журнал

Двухкритериальная идентификация кинетических параметров реакции гидроалюминирования олефинов алкилаланами # 12, декабрь 2013 Б01: 10.7463/1213.0645511

Губайдуллин И. М., профессор, д.т.н. Карпенко А. П., Нурисламова Л. Ф., Савелов А. С.

УДК 519.6

Россия, Уфа, Институт нефтехимии и катализа РАН Россия, МГТУ им. Н.Э. Баумана irekmars@mail.ru apkarpenko@mail.ru nurislamova_lf@mail.ru s savelov@mail.ru

Введение

В настоящее время мировая промышленность ежегодно производит десятки тысяч тонн алюминийорганических соединений (АОС), большая часть которых используется в качестве сокатализаторов в полимеризационных процессах для получения полиэтилена, полипропилена, полиизопренового каучука и т.д. В Институте нефтехимии и катализа РАН (ИНК РАН, г. Уфа) ведутся исследования и разработка оригинальных двухкомпонентных нейтральных каталитических систем, состоящих из соединений металла переменной валентности (титана или циркония) и АОС, которые называются металлокомплексными катализаторами [1, 2]. Промышленное использование указанных катализаторов предполагает знание поведения реагирующей системы при изменении внешних условий (температуры, давления, концентрации реагентов и т.п.), что позволяет оптимизировать режим протекания реакции с целью обеспечения его безопасности и максимизации выхода продукта. Другими словами, для управления технологической установкой, в которой происходит

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

Проблема состоит в том, что химические реакции в присутствии металлокомплексных катализаторов представляют собой весьма сложные, быстро протекающие и далеко не полностью изученные процессы. Поставить и провести натурные химические эксперименты в широком интервале изменений условий проведения реакции - очень затратный, а иногда и опасный процесс. Поэтому необходимо, хотя бы качественно, прогнозировать поведение химической реакции во всех возможных условиях. С этой целью на основе результатов проведенных химических опытов необходимо, прежде всего, построить адекватную математическую модель реакции. Для сложных реакций металлокомплексного катализа эта модель включает в себя от 8 до 60 параметров, иногда отличающихся друг от друга на 20 и более порядков [3].

Работа посвящена задаче идентификации кинетических моделей химических реакций указанного класса, которая относится к классу обратных задач химической кинетики. Известно значительное число работ, посвященных этой задаче [4-6]. Большое число таких работ выполнено в ИНК РАН [7-10]. Как правило, в этих работах задача идентификации ставится как однокритериальная задача оптимизации невязки расчетных и экспериментальных данных. Недостаток такой постановки заключается в том, что в процессе идентификации игнорируется важная априорная информация об особенностях кинетики исследуемой химической реакции. Например, если имеются экспериментальные данные для нескольких температур, при которых проводилась реакция, то такая постановка не учитывает зависимость констант скоростей стадий от температуры по уравнению Аррениуса [11]. Если имеется априорная информация об индукционном периоде того или иного из веществ, участвующего в данной

химической реакции [12], то однокритериальная постановка не позволяет принять в расчет эту информацию.

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

Постановка задачи многокритериальной оптимизации (МКО-задачи) фиксирует вектор варьируемых параметров задачи, множество допустимых значений этого вектора, вектор критериальных функций и, возможно, ограничения на значения компонентов этого вектора. Данная информация позволяет, как правило, выделить не одно решение задачи, но множество таких решений (множество Парето). Поэтому полагаем, как это часто делается в современных публикациях в данной области, что решением МКО-задачи является ее множество Парето. Множество Парето занимает в теории многокритериальной оптимизации исключительное место, поскольку, согласно известному принципу Эджворта-Парето, при «разумном» поведении лица, принимающего решения (ЛИР), выбор решения следует производить на множестве Иарето [13].

Известно большое число методов Парето-аппроксимации, то есть методов построения некоторой конечномерной аппроксимации множества Иарето. Обзор таких методов представлен, например, в работе [14]. В данной работе для решения задачи Парето-аппроксимации используем модифицированный метод адаптивных взвешенных сумм (Adaptive Weighted Sum, AWS) [15]. Базовый вариант метода предложили Рю, Ким и Ван (JH. Ryu, S. Kim, H. Wan) в 2009 г. [16]. Метод основан на аддитивной свертке частных критериев оптимальности. В отличие от классического метода суммы взвешенных критериев (WeightedSum, WS) [14], также использующего такую свертку, метод AWS предполагает адаптацию весовых коэффициентов в процессе итераций на основе информации о текущем положении подобласти поиска. Одной из целей работы является апробация

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

Для сокращения затрат на вычисление значений критериальных функций метод ЛЖ8 использует метамодели этих функций. Рассматриваемые обратные задачи химической кинетики описываются системами обыкновенных дифференциальных уравнений, имеющими высокую жесткость. Отсюда вытекает высокая вычислительная сложность этих задач. Тот факт, что метод ЛЖ8 использует метамодели критериальных функций, позволяет значительно сократить вычислительные затраты (без потери точности Парето-аппроксимации).

В первом разделе работы даем постановку двухкритериальных задач идентификации кинетики химических реакций гидроалюминирования олефинов алкилаланами ИЛЮи2 и С1Л1Ви2. Во втором разделе представляем использованное алгоритмическое и программное обеспечение. В третьем разделе приводим результаты вычислительных экспериментов и их обсуждение. В заключении формулируем основные результаты работы и очерчиваем перспективы ее развития.

1. Постановка двухкритериальных задач идентификация кинетических параметров реакции гидроалюминирования олефинов алкилаланами 1.1. Реакция гидроалюминирования олефинов алкилаланом HAlBu2

Математическая модель реакции представляет собой систему обыкновенных дифференциальных уравнений (ОДУ)

^'Х'1! ^^ — + ^2 ^ Х1Х3 ,

dx2| А — 2&1 Х1 2^2 Х2 + ^3 Х1Х3 ^4 Х2 х

Ах3/— - ^3 .Х1Х3 ^4 Х2 Х3, Ах4/ — ^3 Х1Х3 + ^4 Х2 Х3 .

Здесь и далее х - концентрация реагирующего вещества Х{, мол. доли; к -

константа скорости стадий реакции, 1/мин; t - время, мин.; / е [1:4];

Х1=[Ср27гИ2-С1Л1Би2]2 - димер; Х2=[Ср27гИ2-С1Л1Би2] - мономер;

3

(1)

Х3 =HAlBu2; X4=[Cp2ZrH2•HAlBu2• ClAlBu2] - тригидридный комплекс; Ср=С5Н5.

Константы скорости подчиняются закону Аррениуса

к — к0 ехр

ЯТ

Е 1

или 1п к — 1п к0--а —

0 Я Т

(2)

где к0 и Еа (энергия активации, ккал/моль) - константы, подлежащие оценке; Я - газовая постоянная, ккал/(моль*К); Т - температура, К.

Экспериментальные данные представлены значениями концентраций реагирующих веществ при температурах -65 °С, -60 °С, -55 °С. Обозначим натурный эксперимент, проведенный при первой температуре Т65, при второй температуре - Т60, при третьей - Т55. Введем набор констант реакции

К = ( К(1), К(2), К(3)),

где К(1) =( к(1), к«, кз(1), к«), К(2) =( к(2), к22), к32), к42)), К(3) =( к1(3), к23), кз(3), к43)) -

наборы констант скоростей реакции для экспериментов Т65, Т60, Т55 соответственно.

Во введенных обозначениях из уравнений (1) следует совокупность трех систем ОДУ

Ах-^! ^^ — к + кз( Х1Х3,

Ах2!^^ — 2к 2к/( Х/^ + к^зз Х1Х3 ~к4 Х2Х3

Ах3! — кз( Х1Х3 к^^ Х2 Х3, аХ4/ А — кз Х1Х3 + к_4 Х2 Х3.

Экспериментальные данные предоставлены в виде концентраций

"3'

' е[1:3]. (3)

У0')

/и,

Х

1,и

~(') + Х О')

Л1, 4,иу

(4)

где и ' - порядковый номер экспериментальной точки в ' -ом эксперименте, от' - число наборов концентраций (4) для ' -ого эксперимента.

На константы реакций наложены ограничения

k41) > k{21) > k(1\ j e [1:3], (5)

k™ < k(2) < kf\ i e [1:4], (6)

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

Рассматриваем два критерия оптимальности идентификации - критерий MaXimum Squared Error (MXSE) и критерий Method of Least Squares Error (MLSE).

Критерий оптимальности MXSE имеет смысл максимума квадратов невязок концентраций и равен

MXSE (K) = max (yJ - J, j e[1:3], " e[l: mj.

j,"j } 1

Здесь y"J) расчетное значение величины J).

Критерий MLSE отражает максимальную ошибку аппроксимации закона Аррениуса и представляет собой максимум невязок МНК-оценок констант уравнения Аррениуса:

MLSE(K) = max(s\1) )2, i e[1:4], j e [1:3].

i ,j v '

Здесь s\1) = k\1) - k/1) - отклонение расчетного значения k1) константы от ее МНК-оценки 1).

Двухкритериальная задача параметрической идентификации кинетической модели реакции ставится следующим образом: найти такие допустимые константы реакции ki, i e [1: 4], которые минимизируют критерии оптимальности MXSE (K), MLSE (K).

1.2 Реакция гидроалюминирования олефинов алкилаланом ClAlBu2

Математическая модель реакции представляет собой систему ОДУ

Ах1/А — - - ^10, Ах2/А — - 2^4 - + - w1 + + ^15,

Ах3/А — -- w8, Ах4/А — w5 + w8 - w9, Ах5/А — w5 - w6 - w1 + 2w8 - w9 - w11 - w14 - w15, Ах6/А — w9 - w13, аХ7/А — w9 - w11 - w12, Ах8/А — w6 + w1 - w8 + w10,

Ах9/А — -w1 - w2 - w3 - w10 - w12 + w13 + ш Ахю/ А — w2 - w3 + w10 + w12 + w14 - w15, Аххх!А — w1 + w3 - w13, Ах12/А — w2 + w3,

15

Ах13/ А — - w1 - w-

14

Ах14/А — w1 - w2, Ах15/ А — w13,

где приняты следующие обозначения:

— к13 х13 х9 к18 х1хп , — к14 х14 х9 , — ^8 хю х9 ,

w4 — к4 х2 - к1х1, w5 — к9 х2 х3, w6 — к2 х1х5, w1 — к3 х2 х5.

— к12 х8 х3, — к^р^х^х^

10Л4Л5 :

W■

10

к6 х1х9, — кпх1 х5,

(8)

^12 — к15 х1 х9 ^13 — к16 х6 х11 к19 х15 х9 , ^^14 — к^ц х13х5 к20хю , ^^15 — ^5 хюх5 х2х9 .

Здесь Х1= [Cp2ZrH2•aAlBu2]2, Х2= [Cp2ZrH2•aAlBu2], Х3= СН2СНЯ, Х4= Cp2Zra(CH2CH2R), Х5= HAlBu2, Х6= Bu2Al(CH2CH2R), Х7= Cp2ZrHa, Х8= [Cp2ZrH2•HAlBu2•aAlBu2], Х9= aAlBu2, Х10= [Cp2ZrHa•aAlBu2],

Х11= a2AlBu, Х12= C4H8, Х13= Сp2ZrCl2, Х14= Cp2ZraBu,

Х15= aBuAl(CH2CH2R); Я= C5Hll, C6Hlз, ClHl5, C8Hll, Bu= C4H9, Cp= C5H5.

Экспериментальные данные предоставлены в виде, аналогичном (4):

~')

^ п,

х

15, п,-

Х(') + Х(') Л3, 15,и'

' е [1:3], п' е 1: М' \

На константы реакций наложены ограничения

к3 > к4 > к1,

к9 > к12 > 100.

В реакции гидроалюминирования олефинов кинетическая кривая расходования и накопления наблюдаемых веществ имеет ^-образный вид: в начале процесса наблюдается медленное развитие реакции, которое сменяется периодом ускоренного развития реакции. Индукционным периодом вещества в химической реакции называют временной отрезок [0; ], в течение которого происходит медленное изменение концентрации этого вещества. При превышении верхнего предела отрезка [0; ] происходит резкое увеличение или уменьшение скорости изменения концентрации вещества (рисунок 1).

Рисунок 1 - К определению индукционного периода: у' - скорость изменения концентрации рассматриваемого вещества

Как и в предыдущей задаче, рассматриваем два критерия оптимальности идентификации - критерий Sum Squared Error (SSE) и критерий Induction Time Squared Error (ITSE).

Критерий SSE имеет смысл суммы квадратов невязок концентраций:

SSE (K) =z{y(nj - П,)J е[1:3], n, е [l: m, ].

j,n,

Критерий ITSE представляет собой квадрат невязки времени tind :

ITSE (K) = (tmd - L Г Здесь tind и ~nd - расчетное и экспериментальное значения индукционного периода.

Таким образом, двухкритериальную задачу параметрической идентификации реакции гидроалюминирования олефинов алкилаланом ClAlBu2 (5) формулируем следующим образом: найти такие допустимые константы реакций кг, i е [1:20], которые минимизируют критерии оптимальности SSE (K), ITSE (K).

2. Используемые алгоритмы и программное обеспечение 2.1. Постановка задачи многокритериальной оптимизации

Пусть множеством допустимых значений вектора варьируемых

параметров X является ограниченное и замкнутое множество DX с R)X 1. Положим, что критериальная вектор-функция

F (X) = (f1(X), f2(X),..., fF|( X)) со значениями в критериальном пространстве

RF1 определена в области DX. ЛПР стремится минимизировать в этой области каждый из частных критериев оптимальности /(X), f2(X),..., fF|(X), что условно записываем в виде

min F(X) = F(X*) = F*, (9)

X eDX

* *

где векторы X , F - искомое решение задачи многокритериальной оптимизации. Здесь и далее запись вида Щ, где A - вектор или некоторое

счетное множество, означает размерность этого вектора или мощность множества соответственно.

Вектор-функция F(X) выполняет отображение множества DX во множество DF е {F}, которое называется множеством достижимости. Множество и фронт Парето задачи (9) обозначаем D*X, D* соответственно;

DX с DX, Df с Df . Конечномерные аппроксимации этих множеств обозначаем 0X , 0F соответственно и называем архивными.

2.2. Метод адаптивных взвешенных сумм. Основными в методе A WS являются следующие процедуры.

Инициализация метода. Определяем значения свободных параметров метода: S0 - начальный радиус области доверия (trust region radius); р е (0; 1) - коэффициент сужения этой области; Smin - минимальная величина радиуса. Случайным образом в области DX выбираем центральную точку XC.

Определение центральной точки. На итерации (т +1) центральную точку X^+l = X'c отыскиваем среди точек текущей Парето-аппроксимации 0X = 0X (т), построенной на предыдущей итерации т. С этой целью

X F

элементы текущих множеств 0 , 0 сортируем по возрастанию первого частного критерия f1(X) и представляем в виде линейных списков с прежними наименованиями. Расстояние dj архивной точки F® до

ближайших к ней в списке 0F точек определяет формула

dj =11 F0-1 - F0 || + || F0 - Fj+i ||, j е [1: |0|], (10)

где II - евклидова векторная норма. Алгоритм определения центральной точки X'C использует следующее правило.

1) Если 0| > 2, то полагаем X'C = X0 , где

*

j = ai-g

f л

max d; | X0 g XT

^ je[2: |0|-1] j j J

Здесь Хт - множество точек, использованных в качестве центральных на всех предыдущих итерациях [0: т]. Иными словами, за центральную точку принимаем точку, во-первых, наиболее удаленную от других точек

множества 0х в смысле расстояния (10), и, во-вторых, не использованную на предшествующих итерациях.

2) Если 0| = 2, то с равной вероятностью полагаем X' = Х0 или

Хс = Х 0.

3) Если 0| = 1, то принимается Х'с = X0.

Формирование метамоделей. Метамодель тТ+1( X) = т' (X) представляет собой квадратичную аппроксимацию функции /1 (X) в окрестности точки X1',:

т' (X) = /1(X' ) + ^ШС )(X - X' ) + 2(X - X' )Т И (X' )(X - X' ); т2 (X) = /2(X' ) + У/X' ) (X - X' ) + 2(X - X' )Т И2 (X' )(X - X' ).

Здесь V/ (X'), И1 (X') - вектор градиента и матрица Гессе функции / (X) в точке X'; / = 1,2.

Если 0| > 2, то дополнительно строим метамодели

т'р (X) = Я>т[ (X) + ЛР т2 (X), т'ч (X) = (X) + Л т2 (X), а если 0| = 2 или 0| = 1 - метамодель

т'р (X) = А? т' (X) + Лр т2 (X). В 01 > 2 весовые множители (Лр, ЛЯ) = Лр, (Я, Л2) = Лд определяем по

правилу

лр = Ср((-/2(xс) + /2(X®.-1)), (Ш')-/l(xj>-1))), Л = ^((-/2(x®+1) + /М'с )), (/l(X} +1) - /Х(ГС )));

во | 0 | = 2 - по правилу

Лр = ср ((-/2 (X20) + /2 (X®)), (/1 (X,0) - /1 (X!0)));

в = 1 - по правилу Лр = (0,5, 0,5). Константы cp, cq выбираются таким образом, чтобы обеспечить выполнение условий нормировки tf +Äp = tf + =1.

В работе [16] задача формирования квадратичных метамоделей m' (X), m'2 (X) не раскрыта. Остановимся на этой задаче.

Для формирования квадратичных метамоделей m' (X), m'2 (X) используем центральный композиционный план (ЦКП) с центром в точке X'C [17]. Для упрощения вычислений «звездное» плечо в ЦКП выбираем таким образом, чтобы этот план приобрел свойство ортогональности. В качестве ядра ЦКП используем полный факторный эксперимент. Легко видеть, что в общем случае необходимое число испытаний (вычислений значений критериальной функции) для построения одной метамодели равно

n = nk + 2 |X| + n0,

где nk - число испытаний в точках ядра плана, 2 |X| - число «звездных» точек, n0 = 1 - число испытаний в центре плана. Таким образом, в нашем случае (когда |X| = 2) имеем

n = 22 + 2 • 2 +1 = 9.

Серьезной проблемой при использовании ЦКП является генерация эффективной дробной реплики. Используем для этой цели функции Уолша [17].

Решение оптимизационных задач. Данная процедура предполагает решение задач оптимизации

min m' (X) = m' (Х{ ), min m'2 (X) = m'2 (X2 ), (11)

X eD'C X eD'C

где текущую область доверия D'C определяет формула

D'C = {X | X eDX, X - XC <ö*}.

Если J 0 J > 2, то решения Х[ , X2 позволяют отыскать приближенно оптимальные по Парето точки X'p , X'q , принадлежащие области доверия D'C , путем решения оптимизационных задач

mm m'p (X) = m'p (X'p ), mm m'q (X) = m'q (X'q ). (12)

X eD'C X eD'C

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

В процессе итераций текущий радиус области доверия уменьшаем по

правилу 8Т+Х = р 8Т до достижения минимально допустимой его величины Smin. Новое состояние архивного множества 0X получаем путем добавления в него точек X| , X'2 , X'p , X'q и исключения из полученного набора

доминируемых решений. Аналогично, множество 0F формируем путем добавления в него точек F(X( ), F(X2 ), F(X'p ), F(X'q ).

2.3. Используемые модификации метода адаптивных взвешенных сумм

1) Повышение разнообразия множества архивных точек (модификация AWSM1). Наши исследования показали, что в некоторых случаях метод AWS обеспечивает аппроксимацию только части множества Парето, то есть не обеспечивает достаточное разнообразие множества архивных точек 0X [18]. Суть модификации AWSM1 [15] состоит в решении на каждой итерации в крайних точках архивного множества 0X задач

min m' (X) = m' (X' ), min m2 (X) = m2 (X' );

xedl x eD{

min m' (X) = m' (X[ ), min m'2 (X) = m2 (X'2 ).

x gD2 x sD2

Здесь

D; = {X | X e DX, X - X' <^1}, D2 = {X | X e DX, X - X2 <£J.

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

радиус области доверия D'C на каждой итераций, сокращая в результате

разнообразие архивных точек. Идея модификации AWSM2 состоит в смещении центра области доверия «вглубь» области определения, не изменяя при этом ее радиуса [15].

3) Нейросетевая аппроксимация целевых функций (модификация AWSM 3). Идея данной модификации состоит в построении метамоделей критериальных функций с помощью искусственных нейронных сетей [15]. Поскольку на каждой итерации метода AWS аппроксимация критериальной функции производится в пределах «небольшой» области доверия D'C используем радиально-базисные нейронные сети [15].

3. Результаты экспериментов 3.1. Реакция гидроалюминирования олефинов алкилаланом HAlBu2

Фронт Парето для задачи двухкритериальной параметрической идентификации реакции гидроалюминирования олефинов алкилаланом HAlBu2. (п. 1.1) получен с помощью модифицированного AWS-метода AWSM1 и представлен на рисунке 2. Из рисунка следует, что фронт является непрерывным и выпуклым. Метод AWSM1 обеспечил построение плотной и равномерной аппроксимации фронта Парето.

Рисунок 2 - Фронт Парето для обратной задачи химической кинетики реакции гидроалюминирования олефинов алкилаланом НЛ1Ви2 (выделено

решение, выбранное ЛПР)

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

димеризации и взаимодействие с очередной молекулой НЛ1Ви2. В условиях избытка молекул НА1Ви2 наиболее вероятно взаимодействие комплекса Х2 с молекулой НЛ1Ви2, что приводит к образованию экспериментально идентифицированного тригидридного комплекса Х4. Таким образом, энергия активация Е3 должна быть меньше энергии активации Е2. Энергия активации обратной стадии (ей соответствует константа к2) должна быть больше энергии активации прямой стадии (константа к1). Обратимая стадия каталитического процесса определяет общую скорость и эффективность реакции гидроалюминирования, то есть Е2 > Е1, Е2 > Е3, Е2 > Е4. К тому же, константа скорости обратной стадии должна быть больше скорости прямой стадии. Перечисленным условиям удовлетворяет следующий набор констант:

1п к1 = 14,452

6,960 0,002 Т

1п к2 = 60,505

25,804 0,002 Т

2 545 7 688

1п к3 = 5,232--?-; 1п к4 = 19,007--?-; Т е [208; 218], К.

3 0,002 Т 4 0,002 Т

На рисунке 2 выделено решение, соответствующее данному набору. Рисунок 3 иллюстрирует соответствие расчетных и экспериментальных данных. Зависимость констант скорости всех стадий от температуры подчиняется закону Аррениуса (рисунок 4). В интервале рассматриваемых температур полученные результаты удовлетворяют условиям (5), (6).

а) Т=-650С б) Т=-600С в) Т=-550С

Рисунок 3 - Экспериментальные (•) и расчетные (сплошная линия) концентрации вещества Х1: реакция гидроалюминирования олефинов

алкилаланом НА1Ви2

Рисунок 4 - Зависимости констант скоростей от температуры: реакция гидроалюминирования олефинов алкилаланом НА1Ви2; сплошная линия - закон Аррениуса; • - расчетные точки

Результаты исследования показывают, что в реакции с алкилаланом ИЛ1Би12 самой быстрой является стадия перехода комплекса Х2 в тригидридный комплекс Х4, то есть стадия г4 (рисунок 5). Таким образом, на основании полученной кинетической модели рассматриваемой химической реакции подтверждено, что комплексы Х2 -о- Х1 могут взаимодействовать с молекулой НЛ1Би12 с образованием тригидридного комплекса Х3, который ранее считался неактивным.

Рисунок 5 - Скорости стадий реакции в функции температуры: ri - скорость

i-ой стадии

3.2 Реакция гидроалюминирования олефинов алкилаланом ClAlBu2

Фронт Парето для задачи двухкритериальной идентификации реакции гидроалюминирования олефинов алкилаланом ClAlBu2 представлен на рисунке 6. Из рисунка следует, что этот фронт является почти вертикальным, то есть что имеет место слабая зависимость критерия SSE от варьируемых параметров (в рассматриваемых диапазонах их изменения). Задачи с таким фронтом представляют значительную сложность для классического метода WS [15]. Использованный нами метод AWSM1 обеспечил построение плотной и равномерной аппроксимации фронта Парето.

¡TSE

1 1 ■ 1 « 1 1

1

1 1 1 1

О 0,02 0,04 0,06 0,08

SSE

Рисунок 6 - Фронт Парето для обратной задачи гидроалюминирования

олефинов алкилаланом ClAlBu2

Для всех указанных на рисунке точек фронта Парето полученные наборы констант лежат в пределах допустимых значений и одинаково хорошо описывают экспериментальные данные. Для ЛПР все эти наборы равноценны, так что, с формальной точки зрения, задача является однокритериальной. С другой стороны, в реакции гидроалюминирования олефинов алкилаланом QAlBu2 из всех 15 реагирующих веществ наблюдаются лишь вещества Хз, Х15, точнее говоря, известны лишь концентрации этих веществ в некоторые моменты времени. Эти экспериментальные данные одинаково хорошо описывают целые области пространства кинетических параметров. За счет использования априорной информации об индукционном периоде образования продукта реакции удалось сузить эти области. Другими словами, двухкритериальная постановка задачи позволила отбросить заведомо некорректные решения.

Согласие расчетных и экспериментальных данных для рассматриваемой задачи иллюстрирует рисунок 7. Отклонение расчетных данных от экспериментальных в начале реакции объясняется погрешностью экспериментальных данных. Согласно расчетным данным, время

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

0 100 200 300 400

^ мин

Рисунок 7 - Экспериментальные (•) и расчетные (сплошная линия) концентрации вещества Х3: реакция гидроалюминирования олефинов

алкилаланом С1Л1Би12

Заключение

Предшествующие работы, выполненные в ИНК РАН, позволяют сделать следующие выводы.

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

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

реакции, а энергия активации E и предэкспоненциальный множитель ln (k0),

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

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

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

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

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

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

Работа поддержана грантом РФФИ №12-07-00324-а «Структурная и параметрическая идентификация кинетических моделей реакций нейтрального металлокомплексного катализа» и грантом №12-07-31029 «Идентификация механизма реакции гидроалюминирования олефинов параллельными методами».

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

1. Parfenova L.V., Vil'danova R. F., Pechatkina S. V., Khalilov L.M., Dzhemilev U.M. Zr,Al-complexes as new reagents for olefin hydrometallation // J. Organomet. Chem. 2007. Vol. 692. Р. 3424-3429.

2. Pankratyev E.Y., Tyumkina T.V., Parfenova L.V., Khalilov L.M., Khursan S.L., Dzhemilev U.M. DFT study on mechanism of olefin hydroalumination by XAlBui2 in the presence of Cp2ZrCl2 catalyst. I. Simulation of intermediate formation in reaction of HAlBui2 with Cp2ZrCl2 // Organometallics. 2009. Vol. 28, no. 4. P. 968-977.

3. Григорьева Н.Г., Джемилев У.М., Кутепов Б.И., Балаев А.В., Губайдуллин И.М., Хазипова А.Н., Галяутдинова Р.Р. Разработка кинетической модели димеризации а-метилстирола на цеолите типа Y // Химическая промышленность. 2004. № 9. С. 31-36.

4. Быков В.И., Журавлев В.М. Моделирование и оптимизация химико-технологических процессов: учеб. пособие. М.: РХТУ им. Д.И. Менделеева, 2011. 308 с.

5. Ермакова А., Гудков А.В., Аникеев В.И. Идентификация кинетических моделей // Кинетика и катализ. 1997. Т. 38, № 2. С. 309-318.

6. Темкин О.Н. О кинетических моделях многомаршрутных реакций в гомогенном металлокомплексном катализе // Кинетика и катализ. 2012. Т. 53, № 3. С. 326-357.

7. Балаев А.В., Парфенова Л.В., Губайдуллин И.М., Русаков С.В., Спивак С.И., Халилов Л.М., Джемилев У.М. Механизм реакции циклоалюминирования алкенов триэтилалюминием в алюмациклопентаны, катализируемой Cp2ZrCl2 // Доклады АН. 2001. Т. 381, № 3. С. 364-367.

8. Parfenova L.V., Balaev A.V., Gubaidullin I.M., Abzalilova L.R., Pechatkina S.V., Khalilov L.M., Spivak S.I., Dzhemilev U.M. Kinetic Model of Olefins Hydrometallation by HAlBui2 and Al Bui3 in the Presence Cp2ZrCl2 Catalyst // Int. J. Chem. Kinet. 2007. Vol. 39, no. 6. Р. 333-339.

9. Губайдуллин И.М., Рябов В.В., Тихонова М.В. Применение индексного метода глобальной оптимизации при решении обратных задач химической кинетики // Вычислительные методы и программирование. Новые вычислительные технологии. 2011. Т. 12. С. 137-145. Режим доступа: http://num-meth.srcc.msu.ru/ (дата обращения 01.12.2013).

10. Ахметов И.В., Бобренёва Ю.О., Губайдуллин И.М., Новичкова А.В. Математическое моделирование сложных химических реакций в присутствии металлокомплексных катализаторов на основе многоядерных вычислительных систем // Системы управления и информационные технологии. 2013. № 2.1 (52) . С. 111-115.

11. Штиллер В. Уравнение Аррениуса и неравновесная кинетика: пер. с англ. М.: Мир, 2000. 176 с.

12. Хилько А.В., Спивак С.И., Губайдуллин И.М., Парфенова Л.В. О математическом моделировании индукционного периода химических

реакций // Системы управления и информационные технологии. 2008. № 1.2 (31). С. 264-267.

13. Ларичев О.И. Теория и методы принятия решений: учебник для ВУЗов. М.: Университетская книга, Логос, 2006. 392 с.

14. Карпенко А.П., Митина Е.В., Семенихин А.С. Популяционные методы аппроксимации множества Парето в задаче многокритериальной оптимизации. Обзор // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 4. Режим доступа:

http://www.technomag.edu.ru/doc/363023.html (дата обращения 16.12.2013).

15. Карпенко А.П., Савелов А.С. Модифицированный метод адаптивных взвешенных сумм в задаче многокритериальной оптимизации // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2013. № 11. DOI: 10.7463/1113.0632468

16. Ryu J.-H., Kim S., Wan H. Pareto front approximation with adaptive weighted sum method in multiobjective simulation optimization // Proc. of the 2009 Winter Simulation Conference (WSC), 13-16 December 2009, Austin. 2009. P. 623-633. Available at: http://www.informs-sim.org/wsc09papers/060.pdf , accessed 01.11.2013.

17. Асатурян В.И. Теория планирования эксперимента. М.: Радио и связь, 1983. 248 с.

18. Карпенко А.П., Савелов А.С., Семенихин А.С. Адаптивный метод взвешенных сумм в задаче Парето-аппроксимации // Наука и образование. МГТУ им. Н.Э. Баумана. Электрон. журн. 2012. № 6. Режим доступа: http://www.technomag.edu.ru/doc/423283.html (дата обращения 16.12.2013).

scientific periodical of the bauman mstu

SCIENCE and EDUCATION

EL № FS77 - 48211. №0421200025. ISSN 1994-040S

electronic scientific and technical journal

Two-criterion identification of kinetic parameters of olefin hydroalumination by HAlBui2 and ClAlBui2 # 12, December 2013 DOI: 10.7463/1213.0645511

Gubaidullin I.M., Karpenko A.P., Nurislamova L.F., Savelov A.S.

Ufa, Institute of Petrochemistry and Catalysis, Russian Academy of Sciences Bauman Moscow State Technical University, 105005, Moscow, Russian Federation

irekmars@mail.ru apkarpenko@mail.ru nurislamova_lf@mail.ru s_savelov@mail.ru

This work was carried out in the context of investigation and development of original two-component neutral catalyst systems made of mixed-valence metal compounds (titanium and zirconium) and organoaluminum compounds. A problem of identification of kinetic parameters of olefin hydroalumination, which belongs to a class of inverse problems of chemical kinetics, was considered in this paper. As a rule, identification problem is formulated as a single-criterion optimization problem for residuals between calculations and experimental data. The disadvantage of such a formulation consists in the fact that important prior information about kinetic features of a chemical reaction under investigation is ignored during identification. A distinct feature and novelty of this work is a new formulation of the identification problem; the problem was considered as a two-criterion optimization task which allowed to take into account both experimental data on the kinetics of the chemical reaction under investigation and prior information about this reaction. Solution to this problem represents the Pareto set. In order to approximate the Pareto set a modified Adaptive Weighted Sum method was used. One of the purposes of this work is to test the modified AWS method while solving real world two-criterion inverse problems of chemical kinetics. Two-criterion identification problem of kinetics of olefin hydroalumination by halbui2 and clalbui2 was formulated in this work. The algorithms and software which were used are described in this paper. The obtained calculation results were analyzed and discussed.

Publications with keywords: multiobjective optimization, Pareto set, parameter identification, adaptive weighted sum method, inverse problems of chemical kinetics, two-component neutral catalyst system, mixed-valence metal, organoaluminum compounds, olefin hydroalumination

Publications with words: multiobjective optimization, Pareto set, parameter identification, adaptive weighted sum method, inverse problems of chemical kinetics, two-component neutral catalyst system, mixed-valence metal, organoaluminum compounds, olefin hydroalumination

References

1. Parfenova L.V., Vil'danova R. F., Pechatkina S. V., Khalilov L.M., Dzhemilev U.M. Zr,Al-complexes as new reagents for olefin hydrometallation. J. Organomet. Chem., 2007, vol. 692, pp. 3424-3429.

2. Pankratyev E.Y., Tyumkina T.V., Parfenova L.V., Khalilov L.M., Khursan S.L., Dzhemilev U.M. DFT study on mechanism of olefin hydroalumination by XAlBui2 in the presence of Cp2ZrCl2 catalyst. I. Simulation of intermediate formation in reaction of HAlBui2 with Cp2ZrCl2 . Organometallics, 2009, vol. 28, no. 4, pp. 968-977.

3. Grigor'eva N.G., Dzhemilev U.M., Kutepov B.I., Balaev A.V., Gubaydullin I.M., Khazipova A.N., Galyautdinova R.R. Razrabotka kineticheskoy modeli dimerizatsii a-metilstirola na tseolite tipa Y [Development of kinetic model of dimerization of a-methylstyrene on Y-type zeolite]. Khimicheskayapromyshlennost', 2004, no. 9, pp. 31-36.

4. Bykov V.I., Zhuravlev V.M. Modelirovanie i optimizatsiya khimiko-tekhnologicheskikh protsessov [Modeling and optimization of chemical engineering processes]. Moscow, Publ. of Mendeleyev University of Chemical Technology of Russia, 2011. 308 p.

5. Ermakova A., Gudkov A.V., Anikeev V.I. Identifikatsiya kineticheskikh modeley [Identification of kinetic models]. Kinetika i kataliz, 1997, vol. 38, no. 2, pp. 309-318. (English translation: Kinetics and Catalysis, 1997, vol. 38, no. 2, pp. 284-292.).

6. Temkin O.N. O kineticheskikh modelyakh mnogomarshrutnykh reaktsiy v gomogennom metallokompleksnom katalize [Kinetic models of multi-route reactions in homogeneous catalysis with metal complexes (a review)]. Kinetika i kataliz, 2012, vol. 53, no. 3, pp. 326-357. (English translation: Kinetics and Catalysis, 2012, vol. 53, is. 3, pp. 313-343. DOI: 10.1134/S0023158412030123 ).

7. Balaev A.V., Parfenova L.V., Gubaydullin I.M., Rusakov S.V., Spivak S.I., Khalilov L.M., Dzhemilev U.M. Mekhanizm reaktsii tsikloalyuminirovaniya alkenov trietilalyuminiem v alyumatsiklopentany, kataliziruemoy Cp2ZrCl2 [The Mechanism of Cp2ZrCl2-Catalyzed Alkene Cycloalumination with Triethylaluminum to Give Alumacyclopentanes]. Doklady AN, 2001, vol. 381, no. 3, pp. 364-367. (English translation: Doklady Physical Chemistry, 2001, vol. 381, is. 1-3, pp. 279-282. DOI: 10.1023/A:1012900227841 ).

8. Parfenova L.V., Balaev A.V., Gubaidullin I.M., Abzalilova L.R., Pechatkina S.V., Khalilov L.M., Spivak S.I., Dzhemilev U.M. Kinetic Model of Olefins Hydrometallation by HAlBui2 and Al Bui3 in the Presence Cp2ZrCl2 Catalyst. Int. J. Chem. Kinet, 2007, vol. 39, no. 6, pp. 333-339.

9. Gubaydullin I.M., Ryabov V.V., Tikhonova M.V. Primenenie indeksnogo metoda global'noy optimizatsii pri reshenii obratnykh zadach khimicheskoy 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, pp. 137-145. Available at: http://num-meth.srcc.msu.ru/ , accessed 01.12.2013.

10. Akhmetov I.V., Bobreneva Yu.O., Gubaydullin I.M., Novichkova A.V. Matematicheskoe modelirovanie slozhnykh khimicheskikh reaktsiy v prisutstvii metallokompleksnykh katalizatorov na osnove mnogoyadernykh vychislitel'nykh sistem [Mathematical modeling of complex chemical reactions under the influence of metal catalysts based on multi-core computing systems]. Sistemy upravleniya i informatsionnye tekhnologii, 2013, vol. 52, no. 2.1, pp. 111-115.

11. Stiller W. Arrhenius Equation andNon-Equlibrium Kinetics. Teubner, Leipzig, 1989. (Russ. ed.: Stiller W. Uravnenie Arreniusa i neravnovesnaya kinetika. Moscow, Mir, 2000. 176 p.).

12. Khil'ko A.V., Spivak S.I., Gubaydullin I.M., Parfenova L.V. O matematicheskom modelirovanii induktsionnogo perioda khimicheskikh reaktsiy [On the mathematical modeling of induction period of chemical reactions]. Sistemy upravleniya i informatsionnye tekhnologii, 2008, no. 1.2 (31), pp. 264-267.

13. Larichev O.I. Teoriya i metodyprinyatiya resheniy [Theory and methods of decisionmaking]. Moscow, Universitetskaya kniga, Logos, 2006. 392 p.

14. Karpenko A.P., Mitina E.V., Semenikhin A.S. Populyatsionnye metody approksimatsii mnozhestva Pareto v zadache mnogokriterial'noy optimizatsii. Obzor [Review: population methods of Pareto set approximation in multi-objective optimization problem]. Nauka i obrazovanie MGTUim. N.E. Baumana [Science and Education of the Bauman MSTU], 2012, no. 4. Available at: http://www.technomag.edu.ru/doc/363023.html , accessed 16.12.2013.

15. Karpenko A.P., Savelov A.S. Modifitsirovannyy metod adaptivnykh vzveshennykh summ v zadache mnogokriterial'noy optimizatsii [Modified method of adaptive weighted sums in the problem of multi-criterion optimization]. Nauka i obrazovanie MGTU im. N.E. Baumana [Science and Education of the Bauman MSTU], 2013, no. 11. DOI: 10.7463/1113.0632468

16. Ryu J.-H., Kim S., Wan H. Pareto front approximation with adaptive weighted sum method in multiobjective simulation optimization. Proc. of the 2009 Winter Simulation Conference (WSC), 13-16 December 2009, Austin, 2009, pp. 623-633. Available at: http://www.informs-sim.org/wsc09papers/060.pdf , accessed 01.11.2013.

17. Asaturyan V.I. Teoriyaplanirovaniya eksperimenta [Theory of experiment planning]. Moscow, Radio i svyaz', 1983. 248 p.

18. Karpenko A.P., Savelov A.S., Semenikhin A.S. Adaptivnyy metod vzveshennykh summ v zadache Pareto-approksimatsii [Adaptive weighted sum method for solving Pareto-approximation problem]. Nauka i obrazovanie MGTU im. N.E. Baumana [Science and Education of the Bauman MSTU], 2012, no. 6. DOI: 10.7463/0612.0423283