Главная » Программы » Обратная задача равновесия » AUTOEQUIL
Оригинал описания с сайта Лаборатории химической термодинамики МГУ:
Кирьянов Ю.А., Николаева Л.С. AUTOEQUIL, Свидетельство о государственной регистрации программы на ЭВМ, 2008г., №200861226
Разработанный комплекс вычислительных программ AUTOEQUIL ориентирован на математическое моделирование сложных физико-химических равновесных систем как известного, так и неопределенного молекулярного состава при исследовании свойств новых комплексонов и экстрагентов. Расчет химических равновесий проводится по идеальной модели законов действующих масс и материального баланса. Термодинамическое обоснование модели обеспечивает стандартизацию и автоматизацию расчета равновесий. Разработанные эффективные алгоритмы решения “обратной задачи“ обеспечивают ее безаварийное решение, автоматизированный выбор начальных оценок констант равновесий и сходимость к истинным оценкам параметров [1],[2]. Комплекс AUTOEQUIL, работает в нескольких режимах программного меню, задаваемых пользователем в карте ввода, написан на языке программирования FORTRAN, требует 32 килобайта оперативной памяти.
- Евсеев А.М., Николаева Л.С. Математическое моделирование химических равновесий Приложение 3. Кирьянов Ю.А., Николаева Л.С., Евсеев А.М. Автоматизированная система математического моделирования химических равновесий с учетом кинетики баланса масс (AUTUEQUIL) // Издательство Московского университета, 1988, 146-186.
- Евсеев А.М., Николаева Л.С., Кирьянов Ю.А. Журнал физической химии, 1988, том LXII, вып.5, стр.1 153-1175
Приводим инструкцию для подготовки данных к AUTOEQUIL, которая идет вместе с программой. В конце инструкции также приведены примеры, которые отображают возможности программы.
Описание программы AUTOEQUIL
(Вариант включает решение прямой и обратной задачи равновесия, автоматическое построение модели)
(IBM PC)Программа AUTOEQUIL позволяет решать широкий круг задач по расчету химических равновесий в растворах:
- расчет равновесных концентраций всех молекулярных форм в растворе при заданных константах равновесия(устойчивости) - прямая задача равновесия.
- определение констант равновесий по эксперимен- тальным значениям функции отклика методом наименьших квадратов - обратная задача равновесия.
- автоматизированное построение математической модели, наилучшим образом описывающей экспериментальные данные, на заданном наборе реакций. Лучшей полагается модель с минимальным значением взвешенной суммы квадратов отклонений рассчитанных и экспериментальных значений функции отклика.
В программе предусмотрены различные виды отклика и формы ввода экспериментальных данных. Особое внимание уделено обработке результатов pH-метрического титрования.
Управлять прохождением задачи пользователь может с помощью набора переключателей без каких-либо изменений в программе, поэтому для работы с программой нет необходимости знать FORTRAN, a достаточно освоить ввод информации.
Подробно с алгоритмом и методами решения, применяемыми в программе, можно ознакомиться в работах:
- Евсеев А.М. , Николаева Л.С., Кирьянов Ю.А., ЖФХ, 1988, Т.62, N 5, С.1153
- Кирьянов Ю.А, Николаева Л.С, Евсеев А.М Деп. ВИНИТИ, N 2725-В , 1987Г.
- Евсеев А.М., Николаева Л.С. Математическое моделирование химических равновесий, М. МГУ, 1988г.
Подготовка информации для работы с программой AUTOEQUIL
Ввод и вывод информации в программе осуществляется с помощью операторов READ и WRITE.
Ввод информации происходит во всех случаях, когда это не оговорено особо, по формату (5X, 20I3) для целых чисел и по формату (5X, 7Е9.5) для действительных.
Если на одной строке надо разместить целое и действительное числа, то используется формат (5X, I3, F9.5) или (5X, F9.5, I3) в зависимости от последовательности целого и действительного чисел. пять позиций перед первым числом в строке можно использовать для пометок (например, нумерация строк в матрице).
Для работы программы надо приготовить следующую информацию:
1) Заголовок варианта на одной строке(перфокарте). Формат (15А4).
2) IR1, L1, M, N, NM, NMAX - целые числа, задающие размерность задачи:
IR1 - общее число экспериментальных точек(IR1<140)
L1 - число оцениваемых натуральных логарифмов констант устойчивости (параметров равновесий). Допустимо значение L1=0. (L1<14)
М - число независимых компонент физико-химической равновестной системы (М<6)
N - общее число молекулярных форм в стартовой модели (N<37)
NM - число откликов. Если NM=2, то отклики должны быть однотипные (IPH=0 или IPH=2, см.пункт 5). NM=1, 2
NMAX - максимальное число молекулярных форм в расширенной матрице стехиометрических коэффициентов (NMAX<37)
Замечание. Размерность задачи может быть увеличена. Для этого надо изменить размеры массивов в операторах описания.3) IRN, IRK, IRSTP - целые числа:
IRN - номер начальной точки для данного варианта
IRK - номер конечной точки для данного варианта
IRSTP - шаг по экспериментальным точкам
Эти параметры позволяют использовать неполный план эксперимента без иных изменений в вводной информации.4) IKIN, IDIS, ISIMPL, ILG - целые переменные, используемые в качестве переключателей для задания различных режимов прохождения задачи:
IKIN = 0, во всех случаях
IDIS =
0, рассчитываются дисперсии оценок логарифмов констант образования из диаганальн0й матрицы линеаризованной модели для равноточных измерений функции отклика с учетом вычисленной оценки дисперсии воспроизводимости (в этом случае используется схема равноточных измерений с единичными весами)
1, рассчитываются дисперсии оценок логарифмов констант из диагональной матрицы линеаризованной модели с учетом значений дисперсии измерения отклика в точках плана(весовая схема)
ISIMPL = 0, решается обратная задача для стартовой модели 1, решается прямая задача для стартовой модели
ILG = 0, вводятся натуральные логарифмы констант устойчивости (см. 15) 1, вводятся десятичные логарифмы констант.5) IPH, IPRIN, ITITR, IAVTO - переключатели:
IPH
= 0, функция отклика имеет вид суммы концентраций помеченных комплексов, нормированной на значение исходной концентрации одного из компонентов
= 1, функция отклика имеет вид Y = -LOG[H+] (pH) 2, функция отклика имеет вид суммы концентраций помеченных комплексов
IPRIN=
<0, краткая печать результатов по модели
=0, полная печать по модели
=I, где I=1, .., М; кроме полной печати по модели выводится распределение компонента I по всем молекулярным формам в процентах.
ITITR
= 0, ввод экспериментальных данных потенциометрического титрования( см. пп.17-24)
= 1, исходные концентрации компонентов вводятся в явном виде( см. 16)
= -1, данные потенциометрического титрования (см. пп 17-24) вводятся из отдельного файла с расширением .tit
IAVTO
= 0, после рассчета стартовой модели осуществляется перебор свободных комплексов с целью построения модели с минимальной суммой квадратов отклонений рассчитанных и экспериментальных значений функции отклика
= 1, рассчитывается только стартовая модель6) IRG, RANG - целое и действительное числа.
IRG - номер компонента, исходная концентрация которого служит основой для получения первых приближений к оценкам констант равновесия свободных комплексов
RANG - процент от исходной концентрации компонента с номером IRG. Концентрация, состовляющая этот процент от исходной, используется для получения оценок констант7) Е, EPFL - два действительных числа, определяющих точность решения прямой задачи.
Е - точность по приращениям(величина порядка 1.E-4) в натуральных логарифмах равновесных концентраций компонентов при решении системы нелинейных алгебраических уравнений.
EPFL - абсолютная точность по разности между левыми и правыми частями алгебраических уравнений(величина порядка 1.E-11) решение прямой задачи считается полученным, если достигнута точность по одному из критериев.8) ВЕTA - величина критерия для принятия решения о значимости свободного комплекса по изменению суммы квадратов отклонений и о необходимости включения его в модель. при 0
9) IZN, ZNA - целое и действительное числа.они определяют правило исключения варьируемых комплексов из модели по концентрации. Если концентрация некоторого варьируемого комплекса во всех точках меньше процента ZNA от исходных концентраций компонентов IZN, IZN+1, .., M , то этот комплекс исключается из модели и бракуется.
10) NBRAK - целое число, указывающее количество комплексов заведомо исключенных из перебора (бракованных комплексов).
11) IBRAK - целочисленный массив длины NBRAK, содержащий номера бракованных комплексов в произвольном порядке. Если NBRAK=0, то строка, соответствующая массиву IBRAK опускается.
12) INM - целочисленный массив размерности NM*(NMAX+1).
Если отклик не pH (IPH не равно 1), то на NM строках приготовить информацию, описывающую этот отклик(отклики) в формате (5Х, 30I1). В первой строке пробить номер компонента, по которому формируется отклик и одномерный массив INM длины NMAX, состоящий из 0 и 1, расположенных в последовательности, соответствующей расположению комплексов в матрице а (см. 15); причем 1 соответствует комплексу, который (при включении его в модель) входит в отклик, а 0, если нет. Для второго отклика строка формируется аналогично.
Если отклик по pH, то такие строки не нужны.13) IMO1 - целочисленный массив размерности N, определяющий стартовую модель. Он содержит номера комплексов (строк) в матрице А, входящих в стартовую модель.
14) IPI1 - целочисленный массив размерности N, определяющий варьируемые комплексы в стартовой модели. содержит номера варьируемых комплексов в матрице А. Если таких комплексов нет (L1=0), то такая строка не нужна.
15) PI1, A - действительные массивы.
А - двумерный массив размерности NMAX*М, содержащий расширенную компонентную матрицу стехиометрических коэффициентов. Первые M строк соответствуют компонентам, причем если есть компонент (H+), то он стоит на первом месте. На M+1 месте стоит комплекс (OH-) (он обязательно должен присутствовать в матрице, если откликом является pH). Дальнейший порядок комплексов не существенен.
PI1 - одномерный массив размерности NMAX, содержащий натуральные (при ILG=0) или десятичные (при ILG=1) логарифмы констант равновесия (устойчивости).
При работе программы используются константы реакций, которые входят в стартовую модель.
Значения констант из PI1, относящиеся к варьируемым комплексам, рассматриваются как начальные приближения, а остальные как уже известные константы. При автоматизированном построении модели значения констант из PI1 на вводе, относящиеся к бракованным и свободным комплексам, игнорируются, а используются оценки, вычисленные в процессе работы алгоритма.Информацию надо расположить на NMAX строках. На каждой строке надо напечатать натуральный логарифм константы образования и строку матрицы А.
Если ITITR=1, то приготовить информацию в соответствии с пунктом 16.
16) Y, W, B1 - действительные массивы:
Y - одномерный массив размерности IR1. Содержит экспериментальные значения функции отклика в точках плана
W - одномерный массив размерности IR1. Содержит ошибки в экспериментальных значениях функции отклика при весовой схеме и 1 при равноточных измерениях. При переходе от весовой схемы к равноточной можно не изменять массив W, а положить переключатель IDIS=0.
B1 - двумерный массив размерности IR1*M. Содержит матрицу исходных концентраций компонентов в точках плана ( матрица плана B1). Для компонента (H+) вводится не исходная концентрация [H+], а разность [H+]-[OH-]. Это выделяет компонент (H+) среди прочих и, если компонент (H+) присутствует в модели, то он стоит на первой позиции в матрице A (см 15). Из вышесказанного следует, что только для компонента (H+) возможны отрицательные значения в матрице плана (щелочная область).
На каждой строке (всего таких строк IR1) печатается: экспериментальное значение функции отклика, ошибка в отклике и начальные концентрации компонент в точке плана. Если отклика два (NM=2), то далее печатается IR1 строк со значениями функции отклика и ошибками в отклике номер два.Если ITITR=0, to обрабатываются экспериментальные данные потенциометрического титрования и информация вводится в соответствии с пунктами 17-24. Если ITITR=-1, to информация вводится из дополнительного файла с расширением .tit
17) IC - целое число, определяющее количество кривых титрования (столько раз должна быть подготовленна информация в соответствии с пунктами 18-24). IC<11.
18) IRI - целое число, указывающее сколько точек в данной кривой титрования.
19) VN - действительное число. Начальный объем титрования (в мл).
20) RKT, IRKT - действительное и целое числа.
RKT - концентрация титранта. Если титрант щелочь, то число положительное; если кислота, то со знаком "-". Единица измерения - моль/литр
IRKT - номер компонента в модели, который входит и в титрант. Например, если Na входит в модель в качестве компонента под номером 3 и мы титруем щелочью NaOH, to IRKT=3. Если такого компонента нет, то IRKT=0.21) CDP - действительное число. Концентрация сильной одноосновной кислоты, добавленной с целью получения более кислой области. Если такая кислота не добавлена, то CDP=0.
22) PR - действительный массив размерности M-1. Содержит число протонов, внесенных с компонентами при приготовлении раствора. Числа вводятся в порядке, соответствующем очередности компонентов, начиная со второго. Используется для рассчета [H+] − [OH-].
23) CTO - действительный массив размерности М-1. Содержит исходные концентрации всех компонентов, начиная со второго.
24) VT, YT, WT - действительные одномерные массивы размерности IRI, содержащие экспериментальную информацию по данной кривой потенциометрического титрования.
VT - содержит добавленный обЪем титранта в каждой экспериментальной точке (мл).
YT - содержит значения pH в экспериментальных точках.
WT - содержит значения ошибок в pH в экспериментальных точках. Для перехода от весовой схемы к равноточным измерениям достаточно задать IDIS=0 (см. 5), не изменяя значений W.
На каждой строке (всего таких строк IRI) надо напечатать: добавленный объем, pH, ошибки в pH.Замечание. Информация полученная в соответствии с пунктами 17-24 используется для расчета матрицы плана B1, а также массивов Y и W (см. 16).
Замечание. В программе используются полные константы устойчивости, выраженные через компоненты системы. Например, в системе с компонентами ( H+ , L− , Me2+) константы соответствуют следующим выражениям:
для (Me)(H)(L):
K = [(ME)(H)(L)] / ([ME] [H] [L])Для (МЕ)(ОH)2(L):
K = [(ME)(OH)2(L)] [H]² / [ME]Запуск программы AUTOEQUIL на счет
Запуск происходит по директиве: AEQUIL после этого с терминала надо задать имя файла без расширения, в котором записана исходная информация. Результаты работы программы выводятся на терминал и записываются в файл с расширением .rez. Предполагается, что исходная информация хранится в файле с расширением .dat. Если информация по кривым титрования необходимо ввести из отдельного файла (при ITITR=-1), то она располагается в файле с расширеним .tit. Пример: исходная информация заносится в файл pr1.dat, a часть информации при необходимости в файл pr1.tit. Результаты будут записаны в файл pr1.rez
Контрольные примеры
Рассмотрим несколько контрольных примеров. при этом нумерацию пунктов для вводной информаци будем приводить в соответствии с нумерацией пунктов при описании ввода.
Пример 1. Титрование органической кислоты h4l щелочью. Концентрация кислоты 5.0Е-4; концентрация щелочи 0.0556 пусть расширенная компонентная матрица стехиометрических коэффициентов выглядит так:
Номер Молекулярная форма Матрица A Значения констант H+ L4- 1 H+ 1. 0. 0.0 2 L4- 0. 1. 0.0 3 OH- -1. 0. -32.23 4 HL 1. 1. ? 5 H2L 2. 1. ? 6 H3L 3. 1. ? 6 H4L 4. 1. ? В качестве стартовой выберем точно определенную модель, состоящую из реакций 1-3. Поставим задачу найти модель, наилучшим образом описывающую кривую титрования (пункт 24). Тогда исходная информация будет выглядеть следующим образом:
1). Титрование кислоты H4L щелочью
2). IR1=25; L1=0; M=2; N=3; NM=1; NMAX=7
3). IRN=1; IRK=25; IRSTP=1
4). IKIN=0; IDIS=0; ISIML=1
5). IPH=1; IPRIN=0; ITITR=0; IAVTO=0
6). IRG=2; RANG=10.0
7). E=1.E-4; EPFL=1.E-11
8). BETA=0.3
9). IZN=2; ZNA=3.0
10). NBRAK=0
13). IMO1= 1; 2; 3;
15).
PI1 A 0.0 1. 0. 0.0 0. 1. -32.23 -1. 0. 0.0 1. 1. 0.0 2. 1. 0.0 3. 1. 0.0 4. 1. 17). IC=1
18). IRI=25
19). VN=50.
20). RKT=0.0556; IRKT=0;
21). CDP=0.0
22). PR=4;
23). CTO=5.0E-4
24).
N Vt pH W 1 0.0 3.262 0.02 2 0.1 3.370 0.02 3 0.2 3.489 0.02 4 0.3 3.635 0.02 5 0.4 3.785 0.02 6 0.5 4.013 0.02 7 0.6 4.262 0.023 8. 0.7 4.511 0.066 9 0.8 4.842 0.073 10 0.9 5.308 0.086 11 1.0 5.867 0.146 12 1.1 5.940 0.303 13 1.2 6.116 0.377 14 1.3 7.161 0.274 15 1.4 8.200 0.168 16 1.5 8.369 0.124 17 1.6 8.962 0.087 18 1.7 9.398 0.077 19 1.8 9.666 0.074 20 1.9 10.120 0.084 21 2.0 10.326 0.074 22 2.1 10.377 0.067 23 2.2 10.709 0.061 24 2.3 10.774 0.054 25 2.4 10.843 0.051 Вариант 1.1. Достаточно записать эту информацию в файл (например, for001.dat) в нужном формате и запустить задачу на счет. В результате получется такая последовательность моделей:
n Строка A Стартовая модель Модель 2 Модель 3 Модель 4 Мoдель 5 Ошибки 5 1 1. 0. 0.0 0.0 0.0 0.0 0.0 2 0. 1. 0.0 0.0 0.0 0.0 0.0 3 -1 0. -32.23 -32.23 -32.23 -32.23 -32.23 4 1. 1. --- 20.48 20.43 20.45 20.45 0.143 5 2. 1. --- --- 34.74 34.49 34.49 0.197 6 3. 1. --- --- --- 44.57 44.57 0.246 7 4. 1. --- --- --- --- БРАК --- Значение SUM 164.02 38.283 3.988 0.2916 0.2916 SUM - сумма квадратов отклонений экспериментальных и рассчитанных значений функции отклика.
Вариант 1.2. Рассчитаем модель 4 по весовой схеме. Для этого изменим в входной информации следующие числа:
2). L1=3; N=6
4). IDIS=1
5). IAVTO=1
10). NBRAK=1
11). IBRAK= 7;
13). IMO1= 1; 2; 3; 4; 5; 6;
14). IPI1= 4; 5; 6;
15). PI1= 0.0; 0.0; -32.23; 20.45; 34.49; 44.57;В результате решения прямой задачи получится следующий результат:
SUM = 18.733;
ошибки в константах: 0.015; 0.07; 0.066;Вариант 1.3. Решить обратную задачу ( положить ISIMPL=0) по весовой схеме. Получим:
SUM = 17.676
точка минимума: 20.417; 34.698; 44.760;Пример 2. Рассмотрим титрование кислоты H6L с исходной концентрацией C=2.Е-3 в присутствии металла М2+ с исходной концентрацией C=2.E-3. Для получния более кислой области добавлена сильная неорганическая кислота. Исходная концентрация - 4.Е-3. Концентрация щелочи - 0.0975
Пусть расширенная компонентная матрица стехиометрических коэффициентов выглядит так:
Номер Молекулярная форма Матрица А Значения констант H+ L4- M2+ 1 H+ 1. 0. 0. 0.0 2 L(4-) 0. 0. 0. 0.0 3 M2+ 0. 0. 0. 0.0 4 (OH)- -1. 0. 0. -32.23 5 HL 1. 0. 0. 10.01 6 H2L 2. 1. 0. 16.85 7 H3L 3. 1. 0. 21.21 8 H4L 4. 1. 0. 26.26 9 M(OH) -1. 0. 1. -11.04 10 МL 0. 1. 1. 11.00 ? 11 MHL 1. 1. 1. 20.00 ? 12 M2L 0. 1. 2. ? Вариант 2.1. Решение прямой задачи для модели, включающей равновесия 1-11. Значения переменных и массивов на вводе примут такой вид:
1). Модель H - L - M; отклик pH
2). IR1=5; L1=2; M=3; N=11; NM=1; NMAX=12;
3). IRN=1; IRK=5; IRSTP=1;
4). IKIN=0; IDIS=0; ISIMPL=1;
5). IPH=1; IPRIN=0; ITITR=0; IAVTO=1;
6). IRG=2; RANG=10;
7). E=0.0001; EPFL=1.E-11
8). BETA=0.3
9). IZN=3; ZNA=3.0;
10).IBRAK=0;
13).IMO1=1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11;
14).IPI1= 10; 11; 15).
Константы PI1 Матрица А 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 -32.23 -1.0 0.0 0.0 10.01 1.0 1.0 0.0 16.85 2.0 1.0 0.0 21.21 3.0 1.0 0.0 26.26 4.0 1.0 0.0 -11.14 -1.0 0.0 1.0 11.00 0.0 1.0 1.0 20.00 1.0 1.0 1.0 0.0 0.0 1.0 2.0 17). IC=1;
18). IRI=5;
19). VN=25;
20). RKT=0.0975; IRKT=0;
21). CDP= 0.004;
22). PR= 6.0; 0.0;
23). CTO= 0.002; 0.002;
24).
Vt pH W 1 1.0 2.08 1.0 2 2.0 2.32 1.0 3 2.8 2.63 1.0 4 3.2 2.86 1.0 5 3.6 3.40 1.0 Результаты:
SUM=0.007553; константы: 11.0(.833); 20.0(.968)Вариант 2.2. Решение обратной задачи с целью уточнения значений констант 10, 11. Положим ISIMPL=0.
Результаты:
SUM= 0.001415;
Константы: 11.27(0.34); 19.79(0.42)Вариант 2.3. Нахождение модели, описывающей экспериментальные данные, считая константы 1-9 известными, а реакции 10-12 возможными в данной системе.
Положим: L1=0; N=9; ISIMPL=1; IAVTO=0;
IMO1= 1; 2; 3; 4; 5; 6; 7; 8; 9;
IPI1 - отсутствует;
Результаты: SUM=0.001137; 11.57(0.34); 20.17(0.40); дополнительно проверена значимость комплекса 12 и он забракован.Пример 3. Рассмотрим трехкомпонентную систему из примера 2. В качестве откликов возьмем:
1. долю лиганда, связанную с металлом,
2. Долю свободного металла. исходные концентрации компонентов зададим в явном виде, т.е введем матрицу B1.Вариант 3.1. Решим прямую задачу при известных константах 1-11. Кроме равновесных концентраций, получим распределение лиганда в процентах по всем молекулярным формам. Положим:
1). Модель H - L - M; Два отклика
2). IR1=5; L1=2; M=3; N=11; NM=2; NMAX=12;
3). IRN=1; IRK=5; IRSTP=1;
4). IKIN=0; IDIS=0; ISIMPL=1;
5). IPH=0; IPRIN=2; ITITR=1; IAVTO=1;
6). IRG=2; RANG=10;
7). E=0.0001; EPFL=1.E-11
8). BETA=0.3
9). IZN=3; ZNA=3.0;
10).NBRAK=1;
11).IBRAK=12;
12).INM 2000000000111 3001000000000
13).IMO1=1; 2; 3; 4; 5; 6; 7; 8; 9; 10; 11;
14).IPI1= 10; 11;
15).
Константы PI1 Матрица А 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 -32.23 -1.0 0.0 0.0 10.01 1.0 1.0 0.0 16.85 2.0 1.0 0.0 21.21 3.0 1.0 0.0 26.26 4.0 1.0 0.0 -11.14 -1.0 0.0 1.0 11.00 0.0 1.0 1.0 20.00 1.0 1.0 1.0 0.0 0.0 1.0 2.1 16).
Отклики Y ошибки W 0.5 1.0 0.6 1.0 0.7 1.0 0.8 1.0 0.85 1.0 0.4 1.0 0.3 1.0 0.25 1.0 0.2 1.0 0.1 1.0 матрица плана B1 .011635 .001923 .001923 .007593 .001852 .001852 .004568 .001799 .001799 .003121 .001773 .001773 .001713 .001748 .001748 Результаты: SUM= 0.016490; ошибки в константах 10 и 11: 5.10; 0.189;
Вариант 3.2. Получить распределение металла в процентах по всем молекулярным формам для варианта 3.1. Надо положить IPRIN=3.