CC BY
165
Вычислительные технологии
Том 16, № 2, 2011
Алгоритм минимизации энергии Гиббса: Расчет химического равновесия*
Ю.В. Заика
Институт прикладных математических исследований КарНЦ РАН,
Петрозаводск, Россия e-mail: zaika@krc.karelia.ru
Предложен вычислительный алгоритм минимизации энергии Гиббса для определения состава смеси в состоянии химического равновесия. Алгоритм применим и для случая многофазных систем. Основное внимание уделено критическим случаям, когда методы оптимизации, основанные на использовании градиента и гессиана, теряют эффективность.
Ключевые слова: численные методы условной оптимизации, свободная энергия Гиббса, химическое равновесие.
1. Постановка задачи
Проблема определения равновесного состава смеси является одной из традиционных в химической термодинамике. Обстоятельный анализ данной темы и библиография к ней содержатся, например, в [1-3]. Что касается эффективных численных методов, то, по-видимому, здесь отсчет следует вести с классической работы [4] (см. также статьи [5, 6] и литературные ссылки к ним). Несмотря на наличие развитого программного обеспечения, вычислительные проблемы остаются (универсальный алгоритм невозможен). По оценкам экспертов (автор ориентировался на обзор возможностей пакета HSC Chemistry) примерно для 10% задач автоматический режим их решения не удовлетворяет возрастающим требованиям и возникает необходимость постоянного совершенствования алгоритмов. Так, например, для идеальной газовой смеси задача определения равновесного состава при фиксированных температуре и давлении состоит в минимизации по переменным ni энергии Гиббса
k k G = n (G?(T) + RT ln(Pni n-1)) ^ min, n = m,
i=1 i=1
ni
ство молей г-й составляющей смеси; G0 — стандартный химический потенциал; R —
TP ло атмосфер). В более общем случае добавится еще одна операция суммирования (по количеству фаз) и под знаком логарифма появятся множители, называемые коэффициентами активности. Не останавливаясь здесь на терминологии и подробностях (за этим
*Работа выполнена при финансовой поддержке РФФИ (грант № 09-01-00439).
следует обратиться к руководствам по физической химии), отметим лишь, что с математической точки зрения при nj — +0 имеем особенность (ln — —то), что неизбежно сказывается на работоспособности методов оптимизации, основанных на использовании градиента и гессиана, И дело не только в том, что по условиям задачи может потребоваться определение концентраций компонентов менееЧо-18 -10-19 [1]. В процессе итераций на допустимом многограннике при большом объеме промежуточных вычислений и количестве переменных, исчисляемом десятками, трудно контролировать влияние на точность решения задачи появлений "почти нулей".
Цель настоящей работы — предложить алгоритм, который учитывает указанную особенность, К элементарным операциям отнесены решение задачи линейного программирования и поиск минимума выпуклой функции на отрезке. Для них имеется множество эффективных алгоритмов (автор пользовался пакетом Seilab 4,1), Чтобы сосредоточиться на основной идее, рассмотрим классическую постановку задачи [2]. Итак, требуется минимизировать выпуклую однородную функцию
g(n) = g(nb... ,nfc) = i=1 n-i(ci + ln(nn-1)),
где ci — постоянные (paвные G0/RT + ln P); n — количество молей г-го компонента смеси; n = (n1,..., nk)T G = {n = 0 | n ^ 0, 1 ^ i ^ k}; n = n1 + ... + nk = ||n||. Верхний индекс T означает транспонирование, нули линейных пространств обозначаем одним символом, || • || — октаэдрическая норма в Rk (сумма модулей компонент вектора), В дальнейшем считаем q < 0 (см, пример и замечания), иначе ниже вместо максимальной скорости убывания целевой функции следует говорить о минимальной скорости изменения. Материальный баланс выражается ограничением AT n = b, Элементы матрицы A = {ßj}kxm _ неотрицательные целые числа (множество которых обозначим через Z0), rank A = m, нулевые строки отсутствуют, bj > 0 (bj G R+ = {r > 0}), 1 ^ j ^ m < k, В скобках после n — компоненты градиента g, (gradg)i — — то при n — +0, В силу ограничений материального баланса сумма молей ограничена как сверху, так и снизу: 0 < nmin ^ n ^ nmax < +то.
Используя мольные доли xi = n^n, запишем задачу в следующей форме:
k k g(n) = nf (x) = n(cTx + ((x)) — min, ((x) = ^^ n n-1 ln(n n-1) = ^^ xi ln Xj,
i=1 i=1
x G S = {x G Rk| хг ^ 0, X1 + ... + Xk = 1}, ATn = b, A G Zkxm, b G Rm
В пространстве {Rk, || • ||} вектор x = n/n имеет единичную длину и определяет направление, По непрерывности доопределяем xi ln xi = 0 при xi = 0 в силу а ln а — —0, а — +0. Функцию ( (и /) при необходимости можно считать определенной не только на множестве S: ((0) = 0 ((n) = ((tn) = ((x), t > 0 n G R^k, Значения ((n) определя-
ni n
Итак, целевая функция имеет специальную структуру: сумма молей компонентов смеси умножается на функцию, зависящую лишь от распределения мольных долей,
2. Грубая оценка минимума и направления спуска
Экстремумы функции ((x) на S находятся аналитически: max ( = 0 min ( = — ln k. Максимум достигается на базисных векторах ei = (0,..., 1,..., 0)T (смесь вырождается
в компонент). Минимум единственный и достигается на равномерном распределении мольных долей жг = 1/к, Итак, целевая функция имеет двустороннюю оценку:
Допустимое множество D = {n | n ^ 0, ATn = b} компактно (выпуклый многогранник) , минимальное значение g* = min g оценивается решением двух задач линейного программирования: g * G [g_,g+], g_ = minL0, g+ = minL0, n G D, Вектор целевой функции L0 отличается от вектора с равномерной поправкой компонент на — ln k.
Рассмотрим задачу f (x) = cTx + <^(x) — min x G S, Поясним ее смысл. Пусть формально сумма П фиксирована и ограничение n=b
Тогда оптимальное распределение долей xj определяется именно задачей f — min. Приведем аргумент в геометрических терминах. Фиксируем единичный направляющий вектор: n0 G S (n0 = 1). На луче n(t) = tn0 (t ^ 0 интерпретируем как время движения) имеем g(n(t)) = t(cTn° + ^(n0)), Производная по t равна f (n°). Целесообразно выбрать направление n0, вдоль которого функция g убывает (f < 0) наискорейшим образом: f (n0) — min n0 G S, Величина g(n) определяется как значение f на направлении n0 = n/n (распределение мольных долей), умножеиное на время t = n движения из нуля вдоль n0 в точку n, Ограничение ATn = b определяет компромисс между стремлением
gD Фиксируем номер наименьшего q. В общем случае это номер одного из наименьших q, который без ограничения общности считаем равным k. Для упрощения изложения далее подобные оговорки опускаем. Выразив в f (x) переменную xk = 1 — xi — ... — xk_i, получим (с учетом а = 0 ^ а ln а = 0) фупкцию F(y), y = (x1,... ,xk-1)T G [0,1]k_1. На множестве (0,1)k-i функция F строго выпукла, поскольку гессиан положительно определен: = 1/x + 1/xk, F" = 1/xk, 1 ^ i, j ^ k — 1, i = j, Стационарная точка в (0,1)k-i будет единственным минимумом F на [0,1]k_1. Приравнивая производные функции F по xi к пулю, приходим к системе линейных уравнений
(1+exp(c1 — Ck))x1 + x2 +... + xfc_1 = 1, ... , x1 +... + xfc_2 +(1 + exp(cfc_1 — ck))xk_1 = 1.
Вычитаем первое уравнение из остальных. Затем последовательно "идем снизу вверх": xk_1 = x1 exp(c1 — ck_1^, ..,, x2 = x1 exp(c1 — c2). Подставляя полученные выражения в первое уравнение, имеем решение
Затем вычисляем жк = 1 — ж1 — ... — жк-:1, ж0 С (0,1)к, Значение жк получается и подстановкой г = к, С ростом доминирования к-го компонента в смысле ск < с < 0 (ехр(сг — ск) ^ 1 Уг = к) имеем ж0 ~ 1/ехр(сг — ск) и в пределе получаем жг = 0, 1 ^ г ^ к — 1, жк = 1, / = ск, Смесь содержит практически один (сильно доминирующий) компонент. Если, например, имеются два сильно доминирующие компонента (ск = с5 < сг < 0 г / {к, я}), то жк = ~ 1/2 / ~ ск, Другая крайность — отсутствие доминирования: С1 = ... = ск = в Тогда ж0 = 1/к, При этом $(п) = П(в + ^(ж)), минимум ^ реализуется па равномерном распределении. Найденная точка ж0 С (0,1)к строгого минимума / па 5 определяет направление п0 наискорейшего убывання $(п).
Теперь примем во внимание ограничение п = Ь, Двигаясь по экстремальному лучу ¿п° (£ > 0 ж0 > 0), пройдем в общем случае мимо допустимо го множества О,
L0(n) ^ g(n) ^ L0(n), L0(n) = cTn — nlnk, L0(n) = cTn.
1 ^ i ^ k — 1.
Поэтому рассмотрим другие варианты направлений движения. Пусть n1, n2 — решения задач линейного программирования (ЛП) n — min n — max n G D, Следовательно, известен диапазон значений общего количества молей n = ||n||. Для решения задачи n* имеем оценку n* G [nmin,nmax], nmin = n1, nmax = n2, 0 < n1 ^ n2 < Двигаясь вдоль направлений, в точку n1 попадаем за минимальное время t = n, в точку n2 — за максимальное. Обратимся к структуре целевой функции: g(n) = nf (ж), Пусть формально ж — фиксированный векторный параметр, f (ж) < 0, Тогда задача g — min эквивалентна n — max n G D. С другой стороны, в силу Атж = b/n уменьшение п ведет к росту компонент вектора ж. Это может оказаться предпочтительнее, поскольку f (ж) = cTж + c < 0 (когда |q| достаточно велики). Определим точку п минимума g(n) на отрезке [n1,n2] = An2 + (1 — A)n^ A G [0,1], и предварительно фиксируем направление n0 = ж0 = n/||n||. Помимо n1, n2 следует рассмотреть точки n3,..., n6 G D экстремумов оценочных функций L°(n), L°(n), Поскольку <^(n*) G [— ln k, 0], то можно использовать и "промежуточные" функции cTn — Asnlnk, n G D, As G (0,1), Максимумы рассматриваем в силу того, что векторы из максимума в минимум в линейном приближении могут претендовать на направления спуска. Минимум из минимумов g на невырожденных отрезках [nj, nj] определяет направление n° = d° (не обязательно единственное),
Итак, имеем векторы n°, d° (f(n°) ^ f(d°) < 0). На луче tn° целевая функция убывает максимально без учета AT n = b, Второе направление позволяет попасть в D (при соответствующем значении t = ||n||), но в общем случае с меньшей (по абсолютной величине) скоростью убывания. Можно выбирать и компромиссные направления: h = An° + (1 — A)d°, A G [0,1]. Логично также добавить равномерное направление u° = (1/k,..., 1/k)T, которое минимизирует функцию <^(ж): h = A1n° + A2d° + A3u°, где Aj ^ 0, A1 + A2 + A3 = 1, У вектор ob n°, d° могут быть нулевые компоненты с равными номерами. Тогда a priori игнорируются некоторые составляющие смеси. Чтобы h G R+ (hi > 0), достаточно брать A3 > 0, принимая u° как стабилизатор: 0 < A3 < 1,
3. Первое приближение
Фиксируем направление h G SПR+, Целесообразно начинать с предпочтения d° (A2 ~ 1) в выпуклой комбинации n°, d°, v°. Двигаемся по лучу th, t ^ 0 t интерпретируем как время. Обозначим через aj G столбцы матрицы A: (aj, n) = bj, 1 ^ j ^ m. Угловые скобки означают скалярное произведение (ж, У) = ж1У1 + ... + ж^, Н = (■,-)1/2 -длина. Если подставить n = th в ограничения, то из t(aj, h) = bj определяются моменты времени tj = bj (aj, h)-1 > 0 1 ^ j ^ m, при которых достигается материальный баланс по каждому элементу. Сначала воспользуемся методом наименьших квадратов:
m
р2 = m-1^ (t(aj, h) — bj) — min ^ t° = tmin = (Ab, h) |ATh|- > 0. j=1
Значения j bj заданы приближенно, поэтому n = t°h может оказаться приемлемым
р
n
линейное многообразие Л = {z G Rk| ATz = b}:
p = n — A(ATA)-1(ATn — b) = t°H + B, H = h — A(ATA)-1ATh, B = A(ATA)-1b.
Здесь H — ортогональная проекция вектора h на Л0 = {z G Rk| ATz = 0} B — перпендикуляр (см., например, [7]), являющийся решением системы ATz = b с минимальной длиной. Если p G R^ то p G D и можно за первое приближение принять n * = p. В общем случае, не исключая появления в результате вычислений отрицательных компонент pj < 0, дополнительно проектируем на неотрицательный ортант (переопределяем Pj := 0). Получаем вектор p0, хотя при этом, строго говоря, n * = p0 / D, Проектирование последовательно на многообразие Л и ортант следует повторить несколько раз.
При более детальном рассмотрении действуем следующим образом. Проектируем движущуюся точку th, t ^ 0 ортогонально на линейное многообразие Л: p(t) = tH + B, Для уточнения значений t учтем свойства матрицы материального баланса A Сумма столбцов aj есть вектор а+ > 0 с положительными компонентами, поэтому у вектора H есть отрицательные компоненты H (a+,H) = 0, В силv ATB = b имеются поло-
Bj B > 0
p(t) ^ 0 (неравенства вида a»t ^ в) даст отрезок [tbt2], на котором p(t) G D, В качестве первого приближения n * берем точку минимума g(p(t)), t G [t1,t2]. Отметим, что
{ t}
p( t) D
n1 = p0
4. Итерационный процесс
Вначале остановимся на раскрытии неопределенности а ln а (а — +0) численно, В оригинале [4] приводятся шесть знаков после запятой, т, е, точнее задачу решать не потребовалось, Условно считаем, что (в конкретной задаче) мольные доли а = xi менее 10_7 представляют лишь теоретический интерес. Фиксируем соответствующее ис-чезающе малое по смыслу задачи е > 0 (пусть е=10_10), В вычислениях полагаем а G (—е,е) а ln а = 0, При необходимости можно взять отрезок ряда для функции а ln а а G [0, е).
Следующий шаг — анализ ситуации (gradg(n))i = ci + lnxi — — то, xi — +0, Образно говоря, алгоритму градиентного типа хочется шагнуть в этом направлении, но чтобы удержаться в рамках ограничений, приходится умножать большие числа на малые с потерей точности вычислений и непредсказуемыми последствиями для сходимости (в практически важных задачах число переменных доходит до десятков). Если
g(n)
D
min (grad g(n * ),n) = (grad g(n * ),n *), n G D [7], Для текущего s-го приближения n* ~ n * находится решение n = n* задачи линейного программирования (gradg(n* ),n) — min, n G D, Если в точке n * критерий оптимальности не выполняется, то вектор n* — nS указывает направление строгого убывания целевой функции g и приближение n *+1 определяется минимумом g на отрезке [n*, n*], причем g(n*+1) < g(n*),
В рассматриваемой задаче неопределенность в правой части критерия оптимальности раскрывается в силу однородности целевой функции g: (gradg(n),n) = g(n). Ориентируясь на локальную аппроксимацию g(n) = cTn + n^(n) ~ cTn + n^(n*) (^(0) = 0, <^(n) = ^(n/n) = <^(x)), рассмотрим задачу линейного программирования
k
L_(n) = cTn + np(n*) = ^ (ci + ^(n*))ni — min, n G D.
i=1
k
По сравнению с вариантом L+ (n) = (gradg(nj),n) = ^(с + ln)n — min диапа-
¿=1
зон изменения коэффициентов линейной формы сужается с несобственного множества [-то, Cj] до равномерной поправки <^(nj) G [— ln k, 0], Для хорошего первого приближения nj итерации на основе линейной аппроксимации L- (последовательный переход nj ^ nj ^ nj+1) могут привести к решению задачи.
Попытаемся расширить возможности аппроксимации, объединив невырожденность формы L- с экстремальными свойствами L+, Ограничим множество возможных значений коэффициентов линейной формы, поставив барьер неограниченному росту нормы градиента. Для этого нужен масштаб скорости. Логично воспользоваться оценкой L0(n) ^ g(n) ^ L0(n), gradL0 = c — lnk (покомпонентно cj — lnk), gradL0 = c, фиксировать максимальную по абсолютной величине скорость убывания V = ck — ln k (ck = min с < 0) и не позволять недоминирующим компонентам двигаться существенно быстрее (локально, в линейном приближении). Чем меньше зазор между гиперплоскостями (max(L0 — L0) = nmax ln k, n G D), тем обоснованнее такое ограничение скорости. Обратим внимание на следующее обстоятельство. При классической линеаризации (форма L+) в слагаемых долевой функции (cj + ln ж)n текущее приближение используется для фиксации нелинейности, зависящей явно лишь от мольной доли (замена переменной ж на значение xji). Будем придерживаться этой схемы с той лишь разницей, что вследствие возможности ln ж — —то сначала выделим функции ж ln ж (разрешимую особенность), оставляя n свободными переменными после подстаповки ^jj вместо ж^ Формализуем приведенные соображения.
Определим вектор с, состоящий из коэффициентов линейной аппроксимирующей формы Ls(n) = cTn, алгоритмически. Предварительно обнулим массив: с = 0 G Rk, Если ж = nj/nj = ej (вырожденное приближение смеси одной составляющей), то полагаем с = с, поскольку локально g(n) = n(стж + <^(ж)) ~ n(стж + ^^j)) = cTn, При этом Ls(n) = cTn = L-(n) = L0(n) — верхняя оценка g(n), Классическая форма L+ не имеет смысла (ln = — то), Далее уже считаем, что среди мольных долей нет единицы.
Шаг 1, Рассмотрим первое слагаемое в функции g(n), явно выделив разрешимую особенность: n1(c1 + lnж1) = n1c1 + Пж1 lnж1; ж1 = n1/n, n = n1 + ... + nk, Если в каждом слагаемом g(n) заменить ж^а жjj; то в сумме получим линейную форму L-(n) с достаточно узким диапазоном значений коэффициентов (отрезки [cj — ln k,cj]). Задача
[—то, cj]
в аппроксимации L+(n): L+(n) = (gradg(nj),n), (gradg)j = cj + lnжJj G [—то, cj], Если = 0 (< e), то в правой части тождества n1(c1 +lnж1) = n1c1 + Пж1 lnж1 осуществляем подстановку значения жJ1 текущего приближения в особенность ж1 ln ж1 вместо переменной ж^, В сил у жJ1 ln жJl — 0 в форме Ls(n) — cTn коэффициент при n1 полагаем равным с1 = сь Пусть жJ1 > 0 (> e). Если выполняется неравенство
& = V(C1 +ln жJl)-1 ^ 1 ^ (grad g(nj))1 ^ V, V = cfc — ln k< 0,
то присваиваем с1 := c1 +lnжJ1 (как и в L+), Выполнение неравенства означает, что первая компонента градиента по абсолютной величине не превышает порога скорости |V определенного потенциально доминирующим компонентом смеси. Остается рассмотреть случай G (0,1). При = +0 = +0) и = 1 значения коэффициента с1
c1 V
переходов жJ1 — +0 ^ с1 — с1; — 1 — 0 ^ с1 — V,
Проведем формальные преобразования, В тождество
n1(c1 + ln x1) = ?n1(c1 + ln x1) + (1 — с )(n1c1 + nx1 ln x1)
подставим значения с = £2, x1 = x*v В правой части получаем выражение ^1Vn1 + (1 — ^2)(n1c1 + nx* 1 ln x* ^.Параметр с выбран именно для согласования при указанных предельных переходах (для этих целей допустимы и с = v1 > 0), При переменной
^фиксируем множ итель с1 = + (1 — ^2)(с1 + x * 1 ln x * 1), Требование выполнено. Дополнительно из-за множителя n = n1 + ... + nk следует для j = 2,..., k переопределить значения Cj := Cj + (1 — с)x* 1 ln x * 1 = (1 — с)x* 1 ln x * 1 (до шага 1 Cj = 0),
Шаг 2, Переходим к слагаемому n2(c2 + lnx2) = n2c2 + nx2lnx2, Если x*2 = 0, то в силу x* 2 ln x * 2 = 0 к определенному на шаге 1 значению с2 добавляем величину с2: с2 := с2 + с2, Пуст ь x* 2 > 0, Есл и £2 = V/(c2 + ln x* 2) ^ 1 (скорость убывания
C2 c2 +
lnx *2 в соответствии с левой частью равенства. Остается рассмотреть случай £2 G (0,1), В тождество
n2 (C2 + ln x2) = сП2(С2 + ln x2) + (1 — с )(n2C2 + nx2 ln x2)
подставляем значения с = £2, x2 = x *2: £2Vn2 + (1 — £|)(n2c2 + nx*2lnx *2), К значению c2 добавляем Дс2 = £2V + (1 — £|)(c2 + x*2 lnx*2). Кроме того, к с1 и Cj (j = 3,..., k) добавляем число (1 — *2 lnx*2, Отметим, что если на шаге 1 c1 = с1 (как в верхней оценке L0(n)), то отрицательная добавка (1 — £f)x *2lnx *2 "в нужном направлении". Аналогичны выкладки для варианта с = £2+V2, v2 > 0,
Продолжая процесс преобразования слагаемых ni (ci+ln xi) = nici + nxi ln xi последовательно, сформируем вектор с и определим аппроксимирующую форму Ls(n) = cTn, Заметим, что помимо Ls = L0 (когда пайдетея x * i = 1) реализуется и нижняя оценка функции g:x * = (1/k,..., 1/k)T ^ Ls = L_ = L0, Когда в се ^ ^ 1 (принятый ориентир |V| скорости убывания не превышен), то Ls = L+ = (gradg(n*),n). Неограниченный рост |ci| исключен, поскольку особенность а ln а выделена явно и ограничена, Форми-
C grad g
множеетва линейных форм ({L_} С {Ls} С {L+}) достигнута (но предложенный фора ln а
Далее переход от текущего приближения n * ~ n * к следующему стандартен. Решение n*+1 задачи Ls(n) = cTn — min n G D, соединяем отрезком с вектором nj. Приближение n*+1 определяется минимумом функции g(An*+1 + (1 — A)n*), A G [0,1],
{xi }
разно дополнительно рассмотреть отрезки, соединяющие n*+1 с argmax Ls и, например, с n* Остановимся на следующем варианте. Находим минимумы m1,m2 целевой функции g(n) на отрезках [argmin Ls, argmaxLs], [argminLs,n*], Следующее приближение n *+1 выбираем как arg min g(n), n G [m1,m2]. Критерием остановки может служить min Ls(n) = Ls(n*), n G D, или достижение заданной относительной погрешности: (g(n*+1) — g(n*))/g+•_ < 5, (g(n*+1) — g(n *))/(g _ — g+) < 5.
Отметим, что если на некоторой итерации n*+1 = m1 G [argmin Ls, argmaxLs] С D,
D
В любом случае двигаемся из n * в направлении D, поэтому начального включения n * G D нет необходимости придерживаться слишком строго. Если при n * / D оказалось n*+1 = n* , то сдвигаем n* внутрь отрезка [arg min Ls,n*] (в сторону D, arg min Ls G D) и повторяем цикл либо переходим в допустимое множество D: n*+1 = m1.
Квадратичное приближение. В случае большой размерности задачи можно наращивать влияние нелинейности методом продолжения по параметру:
k
g(n, т) = ^^ (cjnj + rnj ln(njn-1^ — min, ATn = b, n G R, т = т1 > 0,...,т = 1.
j=1
Целесообразна также следующая декомпозиция. Заметим, что n о {n, ж}. Это позволяет перейти к формально скалярной экстремальной задаче n^(n) — min n G [nmin, nmax]. Здесь значения функции ^(n) заданы алгоритмически как решения выпуклой задачи сепарабельпого программирования с параметром n:
f (ж) = стж + <^(ж) — min, Атж = bn-1, ж1 + ... + жk = 1, жj ^ 0.
График функции у(а) = а ln а похож на параболу. В силу этого в зависимости от приближения nj слагаемые жj ln ж^ ^ ^такции <^(ж) аппроксимируются параболами: если компонента вектора nj те слишком близка к граничным точкам отрезка [0,1], то парабола строится по трем точкам, иначе используем информацию y(1/e) = — 1/e, y'(1/e) = 0, Промежуточное приближение nj+1 определяется приближенным решением задачи сепарабельного квадратичного программирования с последующими итерациями проектирования на линейное многообразие Л и неотрицательный ортант.
Многофазная задача. В этом случае принципиальных изменений схема вычислений не претерпевает. Целевая функция аддитивна: n = (п(1)т, ..., п(г)т)т, g(n) = g(1)(n(1)) + ... + g(r)(n(r)), Если j-я фаза однокомпонентна, то формально полагаем ж^') = 1 h(j) = 1, скалярная переменная n(j) входит в общую целевую функцию g(n) линейно. Однако подчеркнем, что при формировании распределения мольных долей
n
основная проблема — размерность задачи.
n
применения представленного алгоритма на основе корректировки линейных аппроксимаций. Основной целью этого алгоритма теперь считаем определение общей суммы молей смеси n = n(1) + ... + n(r). Векторная составляющая n(j) текущего приближения позволяет "вырезать" из ограничения Атn = b соответствующую часть Ar(j)n(j) = b(j). Если среди компонент b(j) вектор а b(j) есть нулевые, то понижаем размерность подзадачи, удаляя j-ю строку и столбцы с ненулевыми AjS^- Делим на n(j) (фиксированное число) и "раскрепощаем" мольные доли. Далее перераспределяем мольные доли в рамках каждой многокомпонентной фазы в соответствии с сепарабельпым квадратичным ал-
n
линейной аппроксимации системы в целом и квадратичной аппроксимации "внутри" многокомпонентных фаз).
5. Пример
В качестве примера рассмотрим задачу из [4], решение которой другими методами известно, что позволяет сравнить результаты. Данные представлены в таблице.
Третья составляющая смеси (Н20) является сильно доминирующей: направление n0 = ж0 (ж° + ... + ж! = 1) наискорейшего убывания целевой функции g(n) без учета ограничений Атп = b имеет третью компоненту 0.999 (ограничимся тремя цифрами
Исходные данные задачи, ln P = 3.932
i Компонент GÖ/RT Н, (1ц N, ai2 о, сцз
1 Н -10.021 1 0 0
2 Н2 -21.096 2 0 0
3 н2о -37.986 2 0 1
4 N -9.846 0 1 0
5 n2 -28.653 0 2 0
6 NH -18.918 1 1 0
7 NO -28.032 0 1 1
8 О -14.640 0 0 1
9 02 -30.594 0 0 2
10 он -26.111 1 0 1
ь3 _ 2 1 1
после точки). Минимумы линейных функций n, L0, L0 достигаются в n1 (третья компонента 1, пятая 0,5, пули не упоминаем), максимумы — в n2 (первая компонента 2, четвертая 1, восьмая 1), Получаем предварительные грубые оценки: n* £ [nmin,nmax] = [1.5, 4], g* £ [g-,g+] = [-49.868, -46.414], Минимум g(n) на отрезке [n1,n2] равен -47.41 и достигается на векторе (0,002, 0, 0,989, 0,01, 0,495, 0, 0, 0,01, 0, 0), нормируя который, получаем d0, Выберем направление спуска h = A1n0 + A2d0 + A3u0, отдавая предпочтение d0 (например, А1 = А3 = 1/30). Метод наименьших квадратов дает t0 = 1.53, причем среднеквадратичная невязка достаточно мала: р = 0.025, Проекция t0h на линейное многообразие Л = {z £ | ATz = 6} равна p = (0.02, 0.004, 0.978, 0.02, 0.48, 0.01,
0.008, 0.01, 0, 0.002) и g(p) = -47.434, Уточнение (на отрезке [t1,t2]) приводит к первому приближению n* = (0.03, 0.026, 0.935, 0.027, 0.462, 0.029, 0.02, 0.015, 0.007, 0.017), g(ni) = -47.466. Далее реализуется итерационный процесс. Первая итерация дает уточнение g* ~ -47.66. Получаем установившийся результат с точностью до тысячных: n* « (0.04, 0.148, 0.783, 0.0015, 0.485, 0.002, 0.028, 0.018, 0.0375, 0.097), g* « -47.76.
6. Замечания
1. Задачи малой размерности достаточно хорошо изучены и допускают наглядную интерпретацию в форме графиков и диаграмм. Изложенный алгоритм нацелен на задачи большой размерности, причем на тот встречающийся на практике случай, когда особенность задачи из-за наличия малых мольных долей (в том числе и в промежуточных подзадачах) становится существенной при большом объеме приближенных вычислений. При этом с учетом однородности целевой функции и линейных ограничений следует предварительно нормировать задачу ^заменить на vn) с тем, чтобы, поделив на v = min bj, прийти к такому диапазону молей химичееких элементов (min6 = 1), при котором легче понять, что следует принять за "нуль". Конечно, точность вычислений и результатов следует согласовать с точностью входных данных, но этот вопрос в математической постановке представляется труднообозримым.
2. Пусть существует Cj > 0. Тогда представим целевую функцию в следующей форме:
g(n) = nß + g0(n), g0(n) = n/0(x), f0(x) = + ^(x), C0i < 0, ß > maxCj. Далее ищем направление наискорейшего убывания g0, как описано выше.
3, Помимо направлений спуска h, регулируемыми параметрами являются порог скорости убывания V (или даже пороги V¿ для каждого компонента) и показатели v > 0 скоростей согласования предельных переходов (с = ) при построении аппроксимирующей линейной формы Ls = cTn,
4, Представленный алгоритм не претендует на высокую скорость сходимости. Метод множителей Лагранжа не используется (возможно, это не недостаток). Предложенную схему вычислений можно реализовать как блок генерации начального приближения в алгоритмах более высокого порядка (обычно требующих хорошего начального приближения с положительными компонентами) с целью избежать вырождения вдали от минимума.
Список литературы
[1] Карпов И.К. Физико-химическое моделирование на ЭВМ в геохимии. Новосибирск: Наука, 1981.
[2] Уэйлес С. Фазовые равновесия в химической технологии: Пер. с англ. Ч. 1, 2. М.: Мир, 1989.
[31 Thermodynamic Equilibria and Extrema / A.N. Gorban, B.M. Kaganovich, S.P. Filippov, A.V. Keiko, V.A. Shamansky, I.A. Shirkalin. Springer, 2006.
[4] White W.B., Johnson S.M., Dantzig G.B. Chemical equilibrium in complex mixtures // J. Chem. Phys. 1958. Vol. 28, No. 5. P. 751-755.
[5] Weber С.F. Convergence of the equilibrium code SOLGASMIX // J. Comput. Phys. 1998. Vol. 145. P. 655-670.
[6] Davies R.H., Dinsdale А.Т., Gisby J.A. et al. MTDATA — thermodynamic and phase equilibrium software from the National Physical Laboratory // CALPHAD. 2002. Vol. 26, No. 2. P. 229-271.
[7] Attetkob A.B., Галкин С.В., Зарубин B.C. Методы оптимизации: Учеб. для вузов. М.: Изд-во МГТУ, 2003.
Поступила в редакцию 28 января 2010 г., с доработки — 20 сентября 2010 г.
