CC BY
52
УДК 544.344:519.688
DOI: 10.18698/0236-3933-2021-3-24-45
РАСЧЕТ РАВНОВЕСНОГО СОСТАВА СЛОЖНЫХ ТЕРМОДИНАМИЧЕСКИХ СИСТЕМ С ИСПОЛЬЗОВАНИЕМ ЯЗЫКА Julia И БИБЛИОТЕКИ Ipopt
Г.В. Белов
gbelov@yandex.ru
ОИВТ РАН, Москва, Российская Федерация
МГУ имени М.В. Ломоносова, Москва, Российская Федерация
Аннотация
Рассмотрены возможности использования библиотеки оптимизации Ipopt для расчета фазового и равновесного составов многокомпонентной гетерогенной термодинамической системы. Приведены две функции, предназначенные для расчета равновесного состава и свойств сложных термодинамических систем, написанные на языке программирования Julia. Указанные функции являются ключевыми в программе, интегрированной с базой данных по термодинамическим свойствам индивидуальных веществ ИВТАНТЕРМО и использованной для проведения тестовых расчетов. Как показали проведенные тестовые расчеты, библиотека Ipopt позволяет определять фазовый и химический составы простых и сложных термодинамических систем с достаточно высоким быстродействием. Использование библиотеки JuMP существенно упрощает подготовку исходных данных для библиотеки Ipopt, поэтому приведенные в статье функции очень компактны. Показано, как можно использовать библиотеку Ipopt, когда неизвестно значение температуры термодинамической системы. Предлагаемый в работе подход применим как для анализа равновесий отдельных химических реакций, так и для расчета равновесного состава сложных химически реагирующих систем. Простота предлагаемых функций позволяет легко интегрировать их в прикладные программы, встраивать в другие библиотеки, использовать в сочетании с более сложными моделями (реальный газ, неидеальные растворы, равновесия с ограничениями) и при необходимости модифицировать. Универсальность языка моделирования JuMP позволяет заменить библиотеку Ipopt на другую без существенной модификации текста программы
Ключевые слова
Термодинамическое равновесие, язык программирования Julia, библиотека Ipopt
Поступила 05.08.2020 Принята 05.10.2020 © Автор(ы), 2021
Введение. Под сложной термодинамической системой в настоящей статье подразумевается многокомпонентная гетерогенная термодинамическая система, в которой возможны химические и фазовые превращения. Будем придерживаться следующего общепринятого определения равновесного состояния: это состояние термодинамической системы, характеризующееся при постоянных внешних условиях неизменностью параметров во времени и отсутствием в системе потоков.
Задача расчета равновесного состава многокомпонентных гетерогенных термодинамических систем является достаточно старой. Первые алгоритмы и программы расчета равновесного состава появились еще в конце 40-х — начале 50-х гг. 20 в. [1-3] и были ориентированы на решение задач ракетной техники. Создание таких программ позволило резко сократить дорогое, сложное и опасное экспериментальное исследование топливных композиций. В первых алгоритмах расчета равновесного состава использовался метод констант равновесия, при этом необходимо было решать систему нелинейных уравнений. Позже появились алгоритмы и программы расчета равновесия, в основу которых положены методы минимизации термодинамических потенциалов [4-10]. Поскольку при решении системы нелинейных уравнений зачастую возможно существование нескольких решений, метод констант равновесия не получил широкого распространения. На практике для проведения термодинамических вычислений гораздо чаще применяются методы минимизации энергии Гибб-са или Гельмгольца. С подробным обзором методов и алгоритмов расчета равновесия, созданных до начала 1980-х гг., можно ознакомиться в [11].
Для расчета равновесного состава необходимо располагать информацией о термодинамических свойствах веществ. Поскольку первые программы расчета равновесия были предназначены для анализа продуктов сгорания топливных композиций, возникла необходимость создания справочников по термодинамическим свойствам веществ в широком диапазоне температур. В частности были созданы справочники ТСИВ [12] и JANAF [13]. В дальнейшем информация из этих справочников перенесена в базы данных, с которыми сопряжены наиболее универсальные программные комплексы [14-16]. Успехи в создании и использовании указанных программ, а также создание справочников по термодинамическим свойствам веществ стимулировали внедрение методов химической термодинамики в различные области науки и техники, в частности в металлургию, геохимию, плазмохимию, энергетику. В последнее время можно наблюдать рост популярности расчетов равновесия в области астрофизики, см., например, [17]. Очень популярна разработанная в NASA программа
CEA [18]. В России и странах бывшего СССР широко применяется разработанная Б.Г. Трусовым программа ТЕРРА [19], исходный текст первой версии которой на языке Fortran опубликован в 1981 г. [14].
Трудности решения задачи расчета равновесного состава сложной термодинамической системы обусловлены некоторыми ее особенностями.
Достоверность результатов вычислений существенно зависит от полноты термодинамической базы данных, в частности, база данных должна содержать информацию о возможно большем числе веществ в конденсированном состоянии. При создании модели термодинамической системы информация из базы данных, как правило, загружается автоматически на основании списка реагентов. Однако не все фазы, автоматически включенные в состав модельной термодинамической системы, могут присутствовать в ней в равновесном состоянии, причем заранее не всегда известно, какие это фазы. Поэтому в процессе вычислений необходимо определять не только химический, но и фазовый состав, при этом не должно нарушаться правило фаз.
Вычисленные значения равновесных концентраций веществ варьируются в очень широких пределах, примерно от нескольких десятков и единиц до величин порядка 10-100. Иногда важно знать хотя бы приближенные значения очень небольших равновесных концентраций веществ, например при анализе возможности образования токсичных соединений в продуктах реакций.
Формульная матрица (матрица индексов химических элементов) может быть довольно большой и при этом очень разреженной.
Поэтому, несмотря на давнюю историю создания методов и программ расчета равновесного состава сложных термодинамических систем, время от времени появляются новые разработки [20-25].
В последнее время наблюдаются две тенденции развития средств термодинамического моделирования сложных термодинамических систем.
1. Создание библиотек с открытым исходным кодом для расчета равновесного состава [23, 26], которые можно встраивать в более сложные системы моделирования.
2. Создание средств расчета равновесного состава с использованием дополнительных ограничений, например химико-кинетических [27-31].
Для дальнейшего развития методов и средств термодинамического моделирования целесообразно использовать языки программирования с открытым исходным кодом, например Python или Julia. На наш взгляд, Julia больше подходит для этой цели, поскольку это полноценный язык про-
граммирования, позволяет выполнять расчеты с высоким быстродействием, располагает большим набором библиотек, удобных для решения научных и технических задач, предоставляет удобный интерфейс к библиотекам оптимизации.
В настоящей работе рассмотрена постановка задачи расчета равновесного состава многокомпонентной гетерогенной термодинамической системы, приведен обзор основных подходов к ее решению, рассказано о библиотеках JuMP и Ipopt и приведены тексты двух функций расчета равновесия с использованием этих библиотек. Для удобства изложения ограничимся рассмотрением модели «идеальный газ, идеальный раствор, нулевой объем конденсированных фаз», в соответствии с которой сжимаемость газовой фазы описывается уравнением состояния идеального газа, все растворы являются идеальными, а объем конденсированных веществ пренебрежимо мал.
Использование принципа максимума энтропии. Традиционный подход к решению задачи расчета термодинамического равновесия можно проиллюстрировать на примере использования принципа максимума энтропии, в соответствии с которым энтропия изолированной системы максимальна в состоянии равновесия. Условие изолированности на практике означает, что скорость релаксации процессов внутри системы гораздо больше скорости изменения условий на ее границах.
В данном случае задачу можно сформулировать так:
N
max S (U, V, x), U = const, V = const, X vjixi = bj, j = 1, ..., m,
x eR" i =1 (1)
Xi > 0, i = 1, ..., N,
где U — внутренняя энергия; V — объем; x — неизвестный вектор состава, компонентами которого являются числа молей веществ; N — число веществ в системе; Vji — число атомов химического элемента j в веществе i (формульная матрица); bj — содержание химического элемента j в системе; m — число химических элементов в системе. Если в термодинамическую систему включены газообразные ионы, число химических элементов увеличивается на единицу, суммарный заряд системы предполагается равным нулю.
Для определения координат условного максимума энтропии воспользуемся методом неопределенных множителей Лагранжа [14]. Функция Лагранжа для расчета равновесного состава химически реагирующей системы имеет вид
(
л = s + х hj j=i
N
х VjiXi - bj
i = 1
где X] — неопределенные множители Лагранжа. Неизвестными являются значения А} и хи
Координаты условного максимума энтропии определяются следующей системой уравнений:
— I +Е vji кj = 0, г = 1, •••,N;
Ju, у, xjФi j=i
N
X VjiXi = bj, j = 1, ..., m. i = 1
(2)
(3)
Используя фундаментальное уравнение Гиббса, уравнение (2) можно преобразовать к виду
gi "Z Vj{kj = 0, i = 1,...,N, j = 1
(4)
где gi — энергия Гиббса вещества г, отнесенная к ЯТ (Т — температура, Я — газовая постоянная).
Равенство (4) выполняется для всех веществ, содержание которых в термодинамической системе отлично от нуля. Для остальных веществ выполняется неравенство
gi vjiXj > 0, i = 1,...,N.
j = 1
(5)
Размерность системы нелинейных уравнений (3) и (4) можно существенно сократить, если использовать подстановку:
Xi = yj exp
- gi (T, p) + X Vki X k
k=1
, i e Ij,
(6)
которая следует из (4) для модели «идеальный газ, идеальный раствор». Здесь р — давление; хг — количество вещества г в фазе ]; у} — общее число молей фазы ]; 1] — множество индексов веществ, входящих в раствор ].
Однако в этом случае приходится решать систему нелинейных уравнений, не зная фазового состава системы. Иными словами, решать одновременно две задачи — определение фазового (набора фаз, присутствую-
щих в термодинамической системе) и равновесного составов. Сложность заключается еще и в том, что для произвольно выбранного набора фаз во многих случаях можно рассчитать равновесный состав. При этом получается несколько разных решений, но истинным будет то из них, которое обеспечивает минимальное значение соответствующего термодинамического потенциала. Если система состоит из небольшого числа фаз, задачу можно решить методом перебора, но в общем случае это слишком затратный путь [32].
Можно показать, что неопределенные множители Лагранжа в данной задаче имеют смысл энергии Гиббса химических элементов, например [33], поэтому для равновесного состава выполняется равенство
N т
X хг§1 = Х]. (7)
I=1 }=1
Таким образом, для найденного решения должны выполняться правило фаз Гиббса и соотношения (3)-(7).
Решение системы уравнений (3), (4) можно найти только в том случае, если известны значения температуры и давления или объема. Однако если дополнить систему соотношениями для расчета энтальпии, энтропии, внутренней энергии и уравнением состояния газовой фазы, равновесие можно определить для заданных значений любых двух параметров из списка: температура, давление, объем, энтропия, энтальпия, внутренняя энергия. Именно такой подход реализован в программах ТЕРРА и АСТРА, разработанных Б.Г. Трусовым [14].
Расчет равновесия методом минимизации термодинамического потенциала. К сожалению, приведенный подход нельзя использовать в сочетании с какой-либо универсальной библиотекой оптимизации, поскольку для случая с неизвестной температурой в целевую функцию необходимо включать каким-то образом информацию о температурной зависимости термодинамических функций веществ, включенных в модельную систему, что сделать весьма затруднительно. Гораздо удобнее использовать менее эффективный с вычислительной точки зрения подход, основанный на использовании метода минимизации термодинамического потенциала. Самыми большими плюсами в этом случае являются простота реализации и очень небольшой (несколько строк) размер процедуры расчета равновесного состава, что существенно упрощает возможность использования библиотеки оптимизации в других программах. Главный недостаток — увеличение времени вычислений.
Отметим, что с вычислительной точки зрения при решении задачи расчета равновесия речь идет об определении координат условного экстремума функции нескольких переменных. Число неизвестных при решении прикладных задач может варьироваться от двух до нескольких сотен.
В соответствии с теоремой Дюгема [34] при отсутствии полей равновесное состояние термодинамической системы, исходные массы которой известны, однозначно характеризуется заданными значениями двух термодинамических параметров. Наиболее часто встречающиеся пары параметров: температура-давление (Т, р); температура-объем (Т, V); энтальпия-давление (Н, р) — горение в реакторе проточного типа; энтропия-давление (5, р) — изоэнтропное расширение до заданного давления; внутренняя энергия-объем (и, V) — горение в постоянном объеме; энтропия-объем (5, V) — изоэнтропное расширение до заданного объема.
Формулировка задачи в координатах температура-давление сводится к определению координат условного минимума энергии Гиббса:
N
min G (T, p, x), T = const, p = const, X vjixi = bj, j = 1,..., m,
x e R"
Xi > 0,
i = 1
i = 1,..., N.
(8)
Энергию Гиббса многокомпонентной гетерогенной системы, состоящей из Ыс однокомпонентных конденсированных фаз и Ыш растворов, можно представить в виде
Nc Nm
G (T, p, x) = X GiXi + X i = 1 j = 1
x Xi(Gi + RT ln a)
i e Ij
(9)
где щ — активность вещества г. В безразмерном виде это соотношение можно представить таким образом:
Nc Nm
g(T, p, X) = X giXi + X i = 1 j = 1
x Xi (gi + ln ai)
i e Ij
где ^ = С /(ЯТ), gi = Сг /(ЯТ).
Для модели «идеальный газ, идеальный раствор, нулевой объем конденсированных фаз» справедливо выражение
Nc Nm
g(T, p, X) = X gi Xi +x i = 1 j = 1
X Xi (gi + ln Xi) - yi ln yj
i e Ij
, yj = X Xi; (10)
i eIj
для веществ в конденсированном состоянии —
g? =[H°(T) - TS°(T)]/(RT); для газообразных веществ —
gi = [Hi(T) - TSP(T)]/(RT) + ln(p/ро),
где ро — стандартное давление; Hf(T), Sf(T) — значения энтальпии и энтропии вещества i в стандартном состоянии при температуре T.
Гораздо реже используется формулировка условий равновесия в координатах температура-объем, в соответствии с которой для расчета равновесия необходимо определить координаты условного минимума энергии Гельмгольца:
N
minF(T, У, х), T = const, V=const, X Vjixi = bj, j = 1,..., m,
xeR" i=1 (11)
Xi > 0, i = 1,..., N.
Энергию Гельмгольца многокомпонентной гетерогенной системы, состоящей из Nc однокомпонентных конденсированных фаз и Nm растворов, можно представить в таком виде:
Nc Nm
F (T, V, x) = X GiXi + X I Xi (Gi + RT ln at) - pV. (12)
i = 1 j = 1 i e Ij
Для модели «идеальный газ, идеальный раствор, нулевой объем конденсированных фаз», полагая, что первый раствор — газовая фаза, множество индексов веществ которой обозначены Ig, запишем
Nc
F(T, V, х) = X FiXi + X Xi (Fi + RT ln Xi)
i = 1
i ein
Nm
+ im
j = 2
x Xi (Fi + RT ln Xi) - RT yj ln y}
i e Ii
или в безразмерном виде
f (T, V, X) = X fi°Xi + x Xi (fi + ln Xi) +
(13)
Nc
X
i = 1
i e I„
Nm
+ z j = 2
X Xi {ff +ln Xi)- yjln yj
i e I,
для веществ в конденсированном состоянии —
f0 =[H°(T) - TS?(T)]/(RT); для газообразных веществ —
RT
fi = [H°(T) - TS°(T)]/ (RT) + ln-— - 1.
роУ
Расчет термодинамических функций веществ в стандартном состоянии. Для расчета значений энергий Гиббса и Гельмгольца веществ в стандартных условиях необходимо знать соответствующие значения энтальпии и энтропии. В [12, 35] информация о термодинамических свойствах индивидуальных веществ приведена в виде таблиц и коэффициентов аппроксимирующего полинома ai, с использованием которых температурные зависимости энтропии и изменения энтальпии при данной температуре можно рассчитать по формулам:
S° (T) = a1 + a2 (ln X+1) - a3 / X2 + 2a5X + 3a6X2 + 4a7X3;
Hо (T) - Hо(0) = T (a2 - 2a3 / X2 - a4 / X + a5X + 2a6X2 + 3a7X3 ),
где X = T/10000. Стандартное давление ро равно 1 атм (101 325 Па) [12]. Значение полной энтальпии можно рассчитать по формуле
Hо (T) = A fH о (298,15) + Hо (T) - Hо (0) - [Hо (298,15) -H о (0)],
где A fH о(298,15) — энтальпия образования вещества при 298,15 K.
В базе данных NIST Chemistry Webbook (webbook.nist.gov) приведены коэффициенты аналогичных соотношений для аппроксимации температурной зависимости термодинамических функций:
S°(T) = A ln (t) + Bt + Ct2/2 + Dt 3/3-E / (2t 2) + G;
H °(T) = H о(298,15) + At + Bt2/2 + Ct 3/3 + Dt4/4-E /t + F - H;
t = T/1000; A, B, C, D, E, F, G, H — коэффициенты полинома. Стандартное давление ро = 1 бар.
При использовании полиномов NASA [16] значения термодинамических функций при данной температуре вычисляются так (ai, bi — коэффициенты):
S°(T) / R = -a{T"2 / 2 - a2T-1 + a3 ln T + aAT + + a5T2/2 + aT 3/3 + a7T 4/4 + b^;
Hо (T) / (RT) = -a{T"2 + a2 (ln T) / T + аъ +aAT /2 + a5T2 /3 + + a6T3 /4 + ayT4/5 + MT.
О библиотеках JuMP и Ipopt. Для реализации процедуры вычислений выбран язык программирования Julia [36] и использованы библиотеки JuMP и Ipopt.
Библиотека JuMP — это язык моделирования с открытым исходным кодом [37, 38], который позволяет пользователям формулировать широкий спектр задач оптимизации (линейные, смешанно-целочисленные, квадратичные, коническо-квадратичные, полуопределенные и др.) с использованием высокоуровневого алгебраического синтаксиса.
Первые языки алгебраического моделирования (AML), предназначенные для представления задач линейного программирования и других задач оптимизации в естественной алгебраической форме, созданы в конце 1970-х гг. Сами языки моделирования не решают задач оптимизации, их цель — передать условия задачи в процедуру оптимизации. В качестве примера можно привести GAMS [39] и AMPL [40], два хорошо известных коммерческих AML, разработка которых началась в 1978 г. и 1985 г. На сегодняшний день AMPL, GAMS и другие подобные коммерческие библиотеки характеризуют уровень развития AML. Они широко используются в промышленности и в академических кругах. Данные языки очень эффективны для решения широкого класса проблем, однако имеют ряд недостатков. В частности, они не позволяют реализовать взаимодействие с процедурой оптимизации во время ее работы, которое дало бы возможность контролировать процесс решения и ускорить работу процедуры в том случае, когда необходимо решить ряд однотипных задач. Кроме того, формулировка задачи с использованием этих языков может быть довольно трудоемкой.
Библиотека JuMP дополняет список языков моделирования с открытым исходным кодом, в который входят, в частности, YALMIP [41], CVX [42] и Pyomo [43]. Эти AML встроены в языки программирования общего назначения и удобны для использования, однако низкое быстродействие таких языков, как MATLAB и Python, не позволяет использовать достоинства этих AML в полной мере. Для решения этой проблемы был создан JuMP, который интегрирован с языком программирования Julia.
Альтернативой AML является формулировка задачи на языке низкого уровня. Если задача достаточно большая, зачастую приходится иметь дело с разреженными матрицами, заполнение ячеек которых довольно трудоемко и сопряжено с появлением ошибок. Для этой цели иногда ис-
пользуют специальные процедуры — генераторы матриц. Использование AML позволяет решить проблему заполнения ячеек матриц.
При решении задач нелинейной оптимизации AML по заданным алгебраическим уравнениям генерируют процедуры для точного вычисления производных, которые непосредственно могут вызывать библиотеки оптимизации. При необходимости для вычисления производных вместо AML можно использовать инструменты автоматического дифференцирования.
Технические задачи, которые должен выполнять AML, можно условно разделить на две части: во-первых, загрузить в память проблему, введенную пользователем, а во-вторых, сгенерировать входные данные, требуемые процедурой оптимизации в соответствии с типом проблемы. Обе эти задачи решает библиотека JuMP, которая использует технику автоматического (или алгоритмического) дифференцирования для вычисления производных выражений, введенных пользователем.
Кроме градиентов функций, процедуры оптимизации часто используют матрицы вторых производных. Такого рода матрицы также могут быть вычислены библиотекой JuMP с использованием техники автоматического дифференцирования.
Для решения задачи расчета равновесия в настоящей работе выбрана свободно распространяемая библиотека с открытым исходным кодом Ipopt (https://coin-or.github.io/Ipopt/index.html), предназначенная для определения координат условного экстремума нелинейной функции многих переменных с использованием метода внутренней точки. Описание реализованного в Ipopt алгоритма приведено в [44].
Реализация алгоритма расчета равновесия на языке Julia для (T, p)-и (T, У)-задач. Приведем тексты функций на языке Julia для расчета равновесного состава с использованием библиотек JuMP и Ipopt. В качестве целевой функции используются энергии Гиббса и Гельмгольца, отнесенные к произведению RT.
На вход функций подается следующая информация: число элементов (m), число веществ (к), число растворов (np), число конденсированных фаз, образованных одним веществом (nc), массив безразмерных значений энергий Гиббса или Гельмгольца веществ (g), массив индексов веществ в фазах-растворах (jj), формульная матрица (A), модель для решателя (model). Вещества в списке упорядочены следующим образом. В его начале находятся конденсированные вещества, не входящие в состав растворов, затем следуют газообразные вещества, далее расположены компоненты конденсированных растворов.
Каждая функция возвращает равновесные концентрации веществ x¡, числа молей фаз y и химические потенциалы элементов Xj (sp). В ряде случаев удается лишь приближенно определить фазовый состав и равновесные концентрации веществ, присутствующих в системе в незначительном количестве. Для уточнения фазового состава и значений равновесных концентраций веществ используются соотношения (5) и (6).
function calc_Gibbs(m,k,np,nc,g,jj,A,b,model)
@variable(model, x[1:k] >= 0, start = 1.e-3) @variable(model, y[1:np] >= 0, start = 1.e-3)
@NLobjective(model, Min, sum(x[i]*g[i] for i in 1:nc)+sum(sum(x[i]*(log(x[i]) + + g[i]) for i in jj[1,j]:jj[2,j]) - y[j]*log(y[j]) for j in 1:np)) for j in 1:np
@constraint(model, sum(x[i] for i in jj[1,j]:jj[2,j]) == y[j]) end
@constraint(model, con, A'*x .== b)
JuMP.optimize!(model)
equilibrium_concentrations=zeros(k)
for i in 1:k equilibrium_concentrations[i]=value(x[i]) end
phase_moles=zeros(np)
for i in 1:np phase_moles[i]=value(y[i]) end
sp=zeros(m)
for i in 1:m sp[i]=shadow_price(con[i]) end return equilibrium_concentrations, phase_moles, sp end
function calc_Helmholtz(m,k,np,nc,g,jj,A,b,model)
@variable(model, x[1:k] >= 0, start = 1.e-3) @variable(model, y[1:np] >= 0, start = 1.e-3) @NLobjective(model, Min, sum(x[i]*g[i] for i in 1:nc)+ sum(x[i]*(log(x[i]) + g[i]) for i in jj[1,1]:jj[2,1])+
sum(sum(x[i] *(log(x[i]) + g[i]) for i in jj[1,j]:jj[2,j]) - y[j]*log(y[j]) for j in 2:np)) for j in 2:np
@constraint(model, sum(x[i] for i in jj[1,j]:jj[2,j]) == y[j]) end
@constraint(model, con, A'*x .== b)
JuMP.optimize!(model)
@show objective_value(model)
equilibrium_concentrations=zeros(k)
for i in 1:k equilibrium_concentrations[i]=value(x[i]) end
phase_moles=zeros(np)
for i in 1:np phase_moles[i]=value(y[i]) end sp=zeros(m)
for i in 1:m sp[i]=shadow_price(con[i]) end return equilibrium_concentrations, phase_moles, sp end
Вызов функции можно осуществить так:
model = Model(Ipopt.Optimizer)
x, y, sp = calc_Gibbs(m,k,np,nc,g,jj,A,b,model)
Реализация алгоритма для случая, когда температура не задана. Если температура не задана, энергии Гиббса и Гельмгольца рассчитать невозможно. В этом случае приходится определять корень T (температуру) нелинейного уравнения:
N
Z - X Xi (T) Zi = 0,
i = 1
где Z — значение заданного параметра (энтальпия, энтропия, внутренняя энергия); x\(T) — равновесные концентрации веществ, отвечающие текущему значению температуры; Zi — парциальное мольное значение параметра (энтальпия, энтропия, внутренняя энергия). Если давление известно, расчет состава проводится с использованием процедуры calc_Gibbs, если задан объем — с использованием процедуры calc_Helmholtz. Для определения корня (температуры) используется библиотека Roots.jl (https://github.com/JuliaMath/Roots.jl).
В качестве примера приведем вызов функции для расчета состава при заданных значениях давления и энтальпии:
T = find1root(find_enthalpy, tmin, tmax, eps),
tmin, tmax — верхняя и нижняя границы интервала для поиска корня; eps — допустимая погрешность расчета; find_enthalpy — функция, в которой по найденным значениям температуры и химического состава вычисляется энтальпия термодинамической системы с использованием информации о термодинамических свойствах веществ. Функция find1root выглядит так:
function find1root(f,a,b,tol) x=find_zero(f,[a,b], atol=tol, Order1()) return x end
Тестирование. Для тестирования работы описанного алгоритма расчета равновесного состава написана программа на языке программирования
Julia (версия 1.5), сопряженная с базой данных ИВТАНТЕРМО [35], которая содержит сведения о термодинамических свойствах более 3300 веществ. Результаты расчетов сравнивались с результатами другой программы расчета равновесного состава [21, 22].
Расчеты выполнены для нескольких относительно простых и более сложных термодинамических систем, в частности для гомогенных газофазных систем с возможностью ионизации и без нее, для гетерогенных систем с газовой фазой и однокомпонентными конденсированными фазами, а также для сложных гетерогенных термодинамических систем с конденсированными растворами и однокомпонентными конденсированными фазами с газовой фазой и без нее. Все задачи решены верно, однако в некоторых случаях решение найти не удалось. Наиболее сложная задача включала в себя список веществ, образованных 22 химическими элементами, при этом из базы данных было отобрано и включено в систему ~ 130 конденсированных веществ, из которых образованы два раствора и несколько десятков однокомпонентных фаз.
Кроме того, выполнены тестовые расчеты c традиционно сложными для численных алгоритмов условиями, при которых ранг матрицы ограничений материального баланса становится меньше числа химических элементов. В качестве примера можно привести системы с одним реагентом, такие как H2O или CO2 в области относительно низких температур, в которых в равновесии присутствует фактически только одно химической соединение, присутствие остальных продуктов реакции пренебрежимо мало. С этими задачами программа также справилась успешно.
Время подготовки программы к расчету, которое включает время, необходимое на загрузку библиотек и базы данных, а также время компиляции программы, занимает примерно минуту (Intel Core Í7-9750H, 2.6 GHz, оперативная память 8 ГБ). При заданной температуре непосредственное время вычислений слабо зависит от размерности задачи и составляет от нескольких сотых до десятых долей секунды. Если температура не задана, время расчета может возрасти до нескольких секунд.
Основные результаты и выводы. Приведены две функции, предназначенные для расчета равновесного состава и свойств многокомпонентных гетерогенных термодинамических систем, написанные на языке программирования Julia. Указанные функции входят в состав программы, сопряженной с базой данных по термодинамическим свойствам индивидуальных веществ ИВТАНТЕРМО, но могут быть включены в состав любой другой аналогичной программы.
Проведенные тестовые вычисления показали хорошую работоспособность функций и сравнительно высокое быстродействие вычислений.
Главными достоинствами предложенного кода являются его открытость, компактность и простота, что позволяет встраивать его в другие приложения, использовать в сочетании с более сложными моделями (реальный газ, неидеальные растворы, равновесия с ограничениями) и при необходимости модифицировать.
Расчет равновесного состава выполнен с использованием библиотеки Ipopt, однако универсальность языка моделирования JuMP позволяет заменить эту библиотеку на другую без существенной модификации текста программы.
ЛИТЕРАТУРА
[1] Brinkley Jr. S.R. Calculation of the equilibrium composition of systems of many constituents. J. Chem. Phys., 1947, vol. 15, no. 2, pp. 107-110.
DOI: https://doi.org/10.1063/L1746420
[2] Kandiner H.J., Brinkley Jr. S.R. Calculation of complex equilibrium relations. Ind. Eng. Chem, 1950, vol. 42, no. 5, pp. 850-855. DOI: https://doi.org/10.1021/ie50485a030
[3] Huff V.N., Gordon S., Morrell V.E. General method and thermodynamic tables for computation of equilibrium composition and temperature of chemical reactions. Report NACA-TR-1037. NASA, 1951.
[4] White W.B., Johnson S.M., Dantzig G.B. Chemical equilibrium in complex mixtures. J. Chem. Phys., 1958, vol. 28, no. 5, pp. 751-755. DOI: https://doi.org/10.1063/L1744264
[5] Синярев Г.Б. Полные термодинамические функции и использование их при расчете равновесных состояний сложных термодинамических систем. Известия вузов. Транспорт и энергетическое машиностроение, 1966, № 2, с. 99-110.
[6] Passy U., Wilde D.J. A geometric programming algorithm for solving chemical equilibrium problems. SIAM J. Appl. Math., 1968, vol. 16, no. 2, pp. 363-373.
[7] Zeleznik F.J., Gordon S. Calculation of complex chemical equilibria. Ind. Eng. Chem, 1968, vol. 60, no. 6, pp. 27-57. DOI: https://doi.org/10.1021/ie50702a006
[8] Алемасов В.Е., Тишин А.П., Дрегалин А.Ф. Общий математический метод расчета и исследования горения при высокой температуре. Физика горения и взрыва, 1971, т. 7, № 1, с. 77-84.
[9] Gordon S., McBride B.J. Computer program for calculation of complex chemical equilibrium compositions, rocket performance, incident and reflected shocks, and Chapman — Jouguet detonations. Report NASA/SP-273. Cleveland, Ohio, Lewis Research Center, 1976.
[10] Karpov I.K., Chudnenko K.V., Kulik D.A., et al. The convex programming minimization of five thermodynamic potentials other than Gibbs energy in geochemical modeling. Am. J. Sc., 2002, vol. 302, no. 4, pp. 281-311.
DOI: https://doi.org/10.2475/ajs.302.4.281
[11] Smith W.R., Missen R.W. Chemical reaction equilibrium analysis: theory and algorithms. John Wiley, 1983.
[12] Гурвич Л.В., Вейц И.В., Медведев В.А. и др. Термодинамические свойства индивидуальных веществ. Т. 1. М., Наука, 1978.
[13] Chase Jr. M.W., Davies C. JANAF thermochemical tables. I, Al-Co. II, Cr-Zr. J. Phys. Chem. Ref. Data, 1985, vol. 14. DOI: https://doi.org/10.18434/T42S31
[14] Синярев Г.Б., Ватолин Н.А., Трусов Б.Г. и др. Применение ЭВМ для термодинамических расчетов металлургических процессов. М., Наука, 1982.
[15] Алемасов В.Е., Дрегалин А.Ф., Тишин А.П. и др. Термодинамические и тепло-физические свойства продуктов сгорания. Т. 1. М., ВИНИТИ, 1971.
[16] McBride B.J., Zehe M.J., Gordon S. NASA Glenn coefficients for calculating thermodynamic properties of individual species. Report NASA/TP-2002-211556. Cleveland, Ohio, Glenn Research Center, 2002.
[17] Agtodez M., Martinez J.I., de Andres P.L., et al. Chemical equilibrium in AGB atmospheres: successes, failures, and prospects for small molecules, clusters, and condensates. Astron. Astrophys, 2020, vol. 637, art. A59.
DOI: https://doi.org/10.1051/0004-6361/202037496
[18] McBride B.J., Gordon S. Computer program for calculation of complex chemical equilibrium compositions and applications. NASA Reference Publication. Cleveland, Ohio, Lewis Research Center, 1994.
[19] Белов Г.В., Трусов Б.Г. Термодинамическое моделирование химически реагирующих систем. М., Изд-во МГТУ им. Н.Э. Баумана, 2013.
[20] Shobu K., Tabaru T. Development of new equilibrium calculation software: CaTCalc. Mater. Trans., 2004, vol. 68, no. 12, pp. 938-987.
DOI: https://doi.org/10.2320/jinstmet.68.983
[21] Белов Г.В. Использование методов линейного программирования для расчета равновесного состава гетерогенных систем с растворами. Вычислительные методы и программирование, 2009, т. 10, № 1, с. 56-61.
[22] Belov G. On linear programming approach for the calculation of chemical equilibrium in complex thermodynamic systems. J. Math. Chem., 2010, vol. 47, no. 1, art. 446. DOI: https://doi.org/10.1007/s10910-009-9580-y
[23] Piro M.H.A., Simunovic S., Besmann T.M., et al. The thermochemistry library Thermochimica. Comput. Mater. Sc., 2013, vol. 67, pp. 266-272.
DOI: https://doi.org/10.1016/j.commatsci.2012.09.011
[24] Ewing M.E., Isaac D.A. Mathematical modeling of multiphase chemical equilibrium. J. Thermophys. Heat Trans., 2015, vol. 29, no. 3, pp. 551-562.
DOI: https://doi.org/10.2514/LT4397
[25] Leal A.M.M., Kulik D.A., Kosakowski G. Computational methods for reactive transport modeling: а Gibbs energy minimization approach for multiphase equilibrium calculations. Adv. Water Resour., 2016, vol. 88, pp. 231-240.
DOI: https://doi.org/10.1016/j.advwatres.2015.11.021
[26] Otis R., Liu Z.K. Pycalphad: Calphad-based computational thermodynamics in Python. J. Open Res. Softw, 2017, vol. 5, no. 1. DOI: http://dx.doi.org/10.5334/jors.140
[27] Keck J.C. Rate-controlled constrained-equilibrium theory of chemical reactions in complex systems. Prog. Energy Combust. Sc., 1990, vol. 16, no. 2, pp. 125-154.
DOI: https://doi.org/10.1016/0360-1285(90)90046-6
[28] Koukkari P., Pajarre R. Calculation of constrained equilibria by Gibbs energy minimization. Calphad, 2006, vol. 30, no. 1, pp. 18-26.
DOI: https://doi.org/10.1016/j.calphad.2005.11.007
[29] Koukkari P. Introduction to constrained Gibbs energy methods in process and materials research. VTT Technical Research Centre of Finland, 2014.
[30] Pajarre R., Koukkari P., Kangas P. Constrained and extended free energy minimization for modelling of processes and materials. Chem. Eng. Sc., 2016, vol. 146, pp. 244258. DOI: https://doi.org/10.1016/jxes.2016.02.033
[31] Ren Z., Lu Z., Gao Y., et al. A kinetics-based method for constraint selection in rate-controlled constrained equilibrium. Combust. Theory Model., 2017, vol. 21, no. 2, pp. 159-182. DOI: https://doi.org/10.1080/13647830.2016.1201596
[32] Белов Г.В. Об определении фазового состава сложных термодинамических систем. ЖТФ, 2019, т. 93, № 6, с. 810-817.
DOI: https://doi.org/10.1134/S0044453719060074
[33] Dorn W.S. Variational principles for chemical equilibrium. J. Chem. Phys., 1960, vol. 32, no. 5, pp. 1490-1492. DOI: https://doi.org/10.1063/L1730947
[34] Пригожин И., Дефэй Р. Химическая термодинамика. Новосибирск, Наука, 1966.
[35] Belov G.V., Dyachkov S.A., Levashov P.R., et al. The IVTANTHERMO — опНш database for thermodynamic properties of individual substances with web interface. J. Phys.: Conf. Ser, 2018, vol. 946, art. 012120.
DOI: https://doi.org/10.1088/1742-6596/946/1/012120
[36] Bezanson J., Edelman A., Karpinsky S., et al. Julia: a fresh approach to numerical computing. SIAM Rev, 2017, vol. 59, no. 1, pp. 65-98.
DOI: https://doi.org/10.1137/141000671
[37] Dunning I., Huchette J., Lubin M. JuMP: a modeling language for mathematical optimization. SIAM Rev., 2017, vol. 59, no. 2, pp. 295-320.
DOI: https://doi.org/10.1137/15M1020575
[38] Legat B., Dowson O., Garcia J.D., et al. MathOptInterface: a data structure for mathematical optimization problems. arxiv.org: веб-сайт.
URL: https://arxiv.org/abs/2002.03447 (дата обращения: 15.04.2021).
[39] Brooke A., Kendrick D., Meeraus A., et al. GAMS: a user's guide. Scientific Press, 1999.
[40] Fourer R., Gay D.M., Kernighan B.W. AMPL: a modeling language for mathematical programming. Thomson/Brooks/Cole, 2003.
[41] Löfberg J. YALMIP: a toolbox for modeling and optimization in MATLAB. Proc. 2004 IEEE Int. Symp. Comp. Aided Control Syst. Des., 2004, pp. 284-289. DOI: https://doi.org/10.1109/CACSD.2004.1393890
[42] Grant M., Boyd S. CVX: MATLAB software for disciplined convex programming. cvXr.com: веб-сайт. URL: http://cvxr.com/cvx (дата обращения: 15.04.2021).
[43] Hart W.E., Laird C., Watson J.P., et al. Pyomo-optimization modeling in Python. Springer, 2017.
[44] Wächter A., Biegler L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., 2006, vol. 106, no. 1, pp. 25-57. DOI: https://doi.org/10.1007/s10107-004-0559-y
Белов Глеб Витальевич — д-р техн. наук, старший научный сотрудник лаборатории теплофизических баз данных ОИВТ РАН (Российская Федерация, 125412, Москва, ул. Ижорская, д. 13, стр. 2); ведущий научный сотрудник лаборатории химической термодинамики кафедры физической химии химического факультета МГУ имени М.В. Ломоносова (Российская Федерация, 119991, Ленинские горы, д. 1, стр. 3).
Просьба ссылаться на эту статью следующим образом:
Белов Г.В. Расчет равновесного состава сложных термодинамических систем с использованием языка Julia и библиотеки Ipopt. Вестник МГТУ им. Н.Э. Баумана. Сер. Приборостроение, 2021, № 3 (136), с. 24-45. DOI: https://doi.org/10.18698/0236-3933-2021-3-24-45
CALCULATION OF EQUILIBRIUM COMPOSITION OF COMPLEX THERMODYNAMIC SYSTEMS USING Julia LANGUAGE AND Ipopt LIBRARY
G.V. Belov gbelov@yandex.ru
Joint Institute for High Temperatures, Russian Academy of Sciences, Moscow, Russian Federation
Lomonosov Moscow State University, Moscow, Russian Federation
Abstract Keywords
The article considers the possibility of using the Ipopt Thermodynamic equilibrium, optimization package for the calculating the phase Julia programming language, and equilibrium compositions of a multicomponent Ipopt package heterogeneous thermodynamic system. Two functions are presented for calculating the equilibrium composition and properties of complex thermodynamic systems, written in the Julia programming language. These functions are the key ones in the program integrated
with the IVTANTERMO database on thermodynamic properties of individual substances and used for conducting test calculations. The test calculations showed that Ipopt package allows determining the phase and chemical compositions of simple and complex thermodynamic systems with a fairly high speed. Using the JuMP modeling language significantly simplifies the preparation of the initial data for the Ipopt package, therefore the functions presented in this article are very compact. It is shown how the Ipopt package can be used when the temperature of the thermodynamic system is unknown. The approach proposed in this work is applicable both for analyzing the equilibrium of individual chemical reactions and for calculating the equilibrium composition of complex chemically reacting systems. The simplicity of the proposed functions allows their easy integrating into application programs, embedding them into more complex applications, using them in combination with more complex models (real gas, nonideal solutions, constrained equilibria), and, if necessary, modifying them. It should be noted that the versatility of the JuMP modeling language makes it possible to replace the Ipopt package with Received 05.08.2020 another one without significant modification of the Accepted 05.10.2020 program text © Author(s), 2021
REFERENCES
[1] Brinkley Jr. S.R. Calculation of the equilibrium composition of systems of many constituents. J. Chem. Phys., 1947, vol. 15, no. 2, pp. 107-110.
DOI: https://doi.org/10.1063/L1746420
[2] Kandiner H.J., Brinkley Jr. S.R. Calculation of complex equilibrium relations. Ind. Eng. Chem, 1950, vol. 42, no. 5, pp. 850-855.
DOI: https://doi.org/10.1021/ie50485a030
[3] Huff V.N., Gordon S., Morrell V.E. General method and thermodynamic tables for computation of equilibrium composition and temperature of chemical reactions. Report NACA-TR-1037. NASA, 1951.
[4] White W.B., Johnson S.M., Dantzig G.B. Chemical equilibrium in complex mixtures. J. Chem. Phys., 1958, vol. 28, no. 5, pp. 751-755. DOI: https://doi.org/10.1063/L1744264
[5] Sinyarev G.B. Full thermodynamics functions and using them at the calculation of complex thermodynamic systems. Izv. vuzov. Transp. i energeticheskoe mashi-nostroenie, 1966, no. 2, pp. 99-110 (in Russ.).
[6] Passy U., Wilde D.J. A geometric programming algorithm for solving chemical equilibrium problems. SIAM J. Appl. Math., 1968, vol. 16, no. 2, pp. 363-373.
[7] Zeleznik F.J., Gordon S. Calculation of complex chemical equilibria. Ind. Eng. Chem., 1968, vol. 60, no. 6, pp. 27-57. DOI: https://doi.org/10.1021/ie50702a006
[8] Alemasov V.E., Tishin A.P., Dregalin A.F. General mathematical method of investigating combustion at high temperatures. Combust. Explos. Shock Waves, 1971, vol. 7, no. 1, pp. 66-71.
[9] Gordon S., McBride B.J. Computer program for calculation of complex chemical equilibrium compositions, rocket performance, incident and reflected shocks, and Chapman — Jouguet detonations. Report NASA/SP-273. Cleveland, Ohio, Lewis Research Center, 1976.
[10] Karpov I.K., Chudnenko K.V., Kulik D.A., et al. The convex programming minimization of five thermodynamic potentials other than Gibbs energy in geochemical modeling. Am. J. Sc., 2002, vol. 302, no. 4, pp. 281-311.
DOI: https://doi.org/10.2475/ajs.302A281
[11] Smith W.R., Missen R.W. Chemical reaction equilibrium analysis: theory and algorithms. John Wiley, 1983.
[12] Gurvich L.V., Veyts I.V., Medvedev V.A., et al. Termodinamicheskie svoystva indi-vidual'nykh veshchestv. T. 1 [Thermodynamic properties of individual materials. Vol. 1]. Moscow, Nauka Publ., 1978.
[13] Chase Jr. M.W., Davies C. JANAF thermochemical tables. I, Al-Co. II, Cr-Zr. J. Phys. Chem. Ref. Data, 1985, vol. 14. DOI: https://doi.org/10.18434/T42S31
[14] Sinyarev G.B., Vatolin N.A., Trusov B.G., et al. Primenenie EVM dlya termo-dinamicheskikh raschetov metallurgicheskikh protsessov [Using computer for thermodynamic calculations of metallurgic processes]. Moscow, Nauka Publ., 1982.
[15] Alemasov V.E., Dregalin A.F., Tishin A.P., et al. Termodinamicheskie i teplofizi-cheskie svoystva produktov sgoraniya. T. 1 [Thermodynamical and thermophysical properties of combustion products. Vol. 1]. Moscow, VINITI Publ., 1971.
[16] McBride B.J., Zehe M.J., Gordon S. NASA Glenn coefficients for calculating thermodynamic properties of individual species. Report NASA/TP-2002-211556. Cleveland, Ohio, Glenn Research Center, 2002.
[17] Agundez M., Martinez J.I., de Andres P.L., et al. Chemical equilibrium in AGB atmospheres: successes, failures, and prospects for small molecules, clusters, and condensates. Astron. Astrophys., 2020, vol. 637, art. A59.
DOI: https://doi.org/10.1051/0004-6361/202037496
[18] McBride B.J., Gordon S. Computer program for calculation of complex chemical equilibrium compositions and applications. NASA Reference Publication. Cleveland, Ohio, Lewis Research Center, 1994.
[19] Belov G.V., Trusov B.G. Termodinamicheskoe modelirovanie khimicheski rea-giruyushchikh system [Thermodynamic modeling of chemically reacting systems]. Moscow, Bauman MSTU Publ., 2013.
[20] Shobu K., Tabaru T. Development of new equilibrium calculation software: CaTCalc. Mater. Trans., 2004, vol. 68, no. 12, pp. 938-987.
DOI: https://doi.org/10.2320/jinstmet.68.983
[21] Belov G.V. Application of linear programming techniques to equilibrium composition calculations for heterogeneous systems with mixtures. Vychislitel 'nye metody i pro-grammirovanie [Num. Meth. Prog.], 2009, vol. 10, no. 1, pp. 56-61 (in Russ.).
[22] Belov G. On linear programming approach for the calculation of chemical equilibrium in complex thermodynamic systems. J. Math. Shem, 2010, vol. 47, no. 1, art. 446. DOI: https://doi.org/10J 007/s10910-009-9580-y
[23] Piro M.H.A., Simunovic S., Besmann T.M., et al. The thermochemistry library Thermochimica. Comput. Mater. Sc., 2013, vol. 67, pp. 266-272.
DOI: https://doi.org/10.1016/j.commatsci.2012.09.011
[24] Ewing M.E., Isaac D.A. Mathematical modeling of multiphase chemical equilibrium. J. Thermophys. Heat Trans., 2015, vol. 29, no. 3, pp. 551-562.
DOI: https://doi.org/10.2514/LT4397
[25] Leal A.M.M., Kulik D.A., Kosakowski G. Computational methods for reactive transport modeling: a Gibbs energy minimization approach for multiphase equilibrium calculations. Adv. Water Resour., 2016, vol. 88, pp. 231-240.
DOI: https://doi.org/10.1016/j.advwatres.2015.11.021
[26] Otis R., Liu Z.K. Pycalphad: Calphad-based computational thermodynamics in Python. J. Open Res. Softw, 2017, vol. 5, no. 1.
DOI: http://dx.doi.org/10.5334/jors.140
[27] Keck J.C. Rate-controlled constrained-equilibrium theory of chemical reactions in complex systems. Prog. Energy Combust. Sc., 1990, vol. 16, no. 2, pp. 125-154. DOI: https://doi.org/10.1016/0360-1285(90)90046-6
[28] Koukkari P., Pajarre R. Calculation of constrained equilibria by Gibbs energy minimization. Calphad, 2006, vol. 30, no. 1, pp. 18-26.
DOI: https://doi.org/10.1016/jj.calphad.2005.11.007
[29] Koukkari P. Introduction to constrained Gibbs energy methods in process and materials research. VTT Technical Research Centre of Finland, 2014.
[30] Pajarre R., Koukkari P., Kangas P. Constrained and extended free energy minimization for modelling of processes and materials. Chem. Eng. Sc., 2016, vol. 146, pp. 244258. DOI: https://doi.org/10.1016/jj.ces.2016.02.033
[31] Ren Z., Lu Z., Gao Y., et al. A kinetics-based method for constraint selection in rate-controlled constrained equilibrium. Combust. Theory Model., 2017, vol. 21, no. 2, pp. 159-182. DOI: https://doi.org/10.1080/13647830.2016.1201596
[32] Belov G.V. Determining the phase composition of complex thermodynamic systems. Russ. J. Phys. Chem, 2019, vol. 93, no. 9, pp. 1017-1023.
DOI: https://doi.org/10.1134/S0036024419060074
[33] Dorn W.S. Variational principles for chemical equilibrium. J. Chem. Phys., 1960, vol. 32, no. 5, pp. 1490-1492. DOI: https://doi.org/10.1063/L1730947
[34] Prigozhin I., Defey R. Khimicheskaya termodinamika [Chemical thermodynamics]. Novosibirsk, Nauka Publ., 1966.
[35] Belov G.V., Dyachkov S.A., Levashov P.R., et al. The IVTANTHERMO — оиИие database for thermodynamic properties of individual substances with web interface. J. Phys.: Conf. Ser, 2018, vol. 946, art. 012120.
DOI: https://doi.org/10.1088/1742-6596/946/17012120
[36] Bezanson J., Edelman A., Karpinsky S., et al. Julia: a fresh approach to numerical computing. SIAM Rev., 2017, vol. 59, no. 1, pp. 65-98.
DOI: https://doi.org/10.1137/141000671
[37] Dunning I., Huchette J., Lubin M. JuMP: a modeling language for mathematical optimization. SIAM Rev., 2017, vol. 59, no. 2, pp. 295-320.
DOI: https://doi.org/10.1137/15M1020575
[38] Legat B., Dowson O., Garcia J.D., et al. MathOptInterface: a data structure for mathematical optimization problems. arxiv.org: website.
Available at: https://arxiv.org/abs/2002.03447 (accessed: 15.04.2021).
[39] Brooke A., Kendrick D., Meeraus A., et al. GAMS: a user's guide. Scientific Press, 1999.
[40] Fourer R., Gay D.M., Kernighan B.W. AMPL: a modeling language for mathematical programming. Thomson/Brooks/Cole, 2003.
[41] Löfberg J. YALMIP: a toolbox for modeling and optimization in MATLAB. Proc. 2004 IEEE Int. Symp. Comp. Aided Control Syst. Des, 2004, pp. 284-289.
DOI: https://doi.org/10.1109/CACSD.2004.1393890
[42] Grant M., Boyd S. CVX: MATLAB software for disciplined convex programming. cvxr.com: website. Available at: http://cvxr.com/cvx (accessed: 15.04.2021).
[43] Hart W.E., Laird C., Watson J.P., et al. Pyomo-optimization modeling in Python. Springer, 2017.
[44] Wächter A., Biegler L.T. On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming. Math. Program., 2006, vol. 106, no. 1, pp. 25-57. DOI: https://doi.org/10.1007/s10107-004-0559-y
Belov G.V. — Dr. Sc. (Eng.), Senior Researcher, Laboratory of Thermophysical Databases, Joint Institute for High Temperatures, Russian Academy of Sciences (Izhorskaya ul. 13, str. 2, Moscow, 125412 Russian Federation); Leading Researcher, Laboratory of Chemical Thermodynamics, Department of Physical Chemistry, Faculty of Chemistry, Lomonosov Moscow State University (Leninskie Gory 1, str. 3, Moscow, 119991 Russian Federation).
Please cite this article in English as:
Belov G.V. Calculation of equilibrium composition of complex thermodynamic systems using Julia language and Ipopt library. Herald of the Bauman Moscow State Technical University, Series Instrument Engineering, 2021, no. 3 (136), pp. 24-45 (in Russ.). DOI: https://doi.org/10.18698/0236-3933-2021-3-24-45
