Оценка точности вычисления минимума функции

[B@R5uk]

ИСН в сообщении #813785 писал(а):

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

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

ИСН, я понимаю, что вы апеллируете к тому, что в природе не существует растворов с двумя устойчивыми состояниями концентраций (хотя, а вдруг!), но одно дело реальный раствор, и совсем другое дело модель.

Кстати, на счёт литературы, alenov, советую вам заглянуть в книжку Сухарёв Тимохов Фёдоров "Курс методов оптимизации". Там много матана и строгих доказательств, но и конкретные алгоритмы тоже приведены (с рекомендациями и с ссылками на другую литературу).

[alenov]

B@R5uk, большое спасибо за подробное решение! Я все это внимательно изучу. С вами удивительно продуктивно общаться!

B@R5uk: "Если же добавить нелинейное ограничение $(x_{OH^-})\cdot(x_{H^+})=\operatorname{const}$ (кстати, чему равна эта константа, если ${x}_{k}$ -- количество вещества в молях?)"

Если константа воды равна 14, то это выражение (ограничение на рН) будет равно 1E-14.
Не спорю ваше решение математически очень изящное и верное. Но оно недостаточно точное для моей цели. Попробуйте добавить последнее условие на рН (в вашей цитате). Фактически это чистая вода и содержание $(x_{H^+})$ должно быть около 1E-7 (рН $= - \log10(x_{H^+})$, в чистой воде при нормальных условиях рН $= 7$). Еще у вас получилось сильно завышенное содержание кислорода. Реально оно не может быть выше 1Е-60. Это следствие: B@R5uk: "Время вычислений с точностью до $10^{-16}$ (максимальная возможная точность)"? Или вы понимаете под максимально возможной точностью отклонение от ограничивающих условий? Если первое, то такая точность абсолютно недостаточна, нужно минимум 1Е-100, если второе ($\sum\limits_{i=1}^N(a_{ij}x_i)=b_j$ - расхождение в этих суммах), то ваш алгоритм работает с достаточной точностью. В настоящий момент я решаю данную задачу за 0.2 секунды (и это время можно еще сократить программными средствами) с суммарным отклонением от ограничивающих условий (1E-8 - значение штрафной функции), это гораздо быстрее, чем в Excel, например. Однако у меня функция Гиббса—Дюгема вычисляется миллионы раз; большая часть ее вычислений производится уже когда решение подобное вашему достигнуто (для этого действительно достаточно порядка 1000 циклов) и нужно добиться требуемой сходимости с ограничивающими условиями.

B@R5uk: "Это условие влияет на концентрации и, следовательно, на количества вещества ионов $H^+$ и $OH^-$. Эти количества вещества связаны через электронейтральность с количеством вещества ионов $HSiO_3^-$, следовательно, если в растворе не будет достаточно много других заряженных ионов, то введение ограничения, связанного с ионным произведением воды, может изменить концентрацию $HSiO_3^-$ на порядок и более."
В литературе по данному вопросу можно найти указание на необходимость условия электронейтральности, про произведение воды (рН) нигде ничего не сказано, я добавил его самостоятельно. По идее, условия электронейтральности может быть достаточно, чтобы привести в порядок баланс ионов в растворе. Именно поэтому я выражаю неуверенность, что оно (условие на рН) необходимо.

B@R5uk: "И всё-таки, откуда взялась функция" Функция Гиббса-Дюгема (Дюхема): http://ru.wikipedia.org/wiki/%D0%A2%D0% ... 0.BC.D0.B0
Для растворов она так записывается, для твердых веществ - потенциал помножается просто на концентрацию. Я к сожалению не могу ответить, почему там можно просто моли вставлять (хотя и должен это знать по роду своего занятия). Можно и концентрации (моль/кг (или литр, что в случае разбавленного раствора непринципиально) - моляльность). Но видимо для самой функции важно только соотношение между компонентами раствора. Попробую найти в специальной литературе.

-- 13.01.2014, 16:25 --

ИСН в сообщении #813785 писал(а):

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

Боюсь, что бывают. При учете неидеальности растворов и введении в эту функцию коэффециентов активности возможность такой ситуации например описана в https://drive.google.com/file/d/0B08REa ... sp=sharing Внизу второго столбца на 590 странице.

B@R5uk, если вы предложите метод, который будет работать быстрее того, который я сейчас использую, например вам удастся расчитать равновесный состав за 10 тыс циклов или меньше при том, что кислород будет меньше 1Е-60 и рН будет около 7 (сейчас 5, что неправдоподобно для чистой воды) - это здорово. Я правда предвижу трудности с воспроизведением алгоритма минимизации, используемого в МатЛабе, потому что в конечном счете сам МатЛаб я использовать не могу. Но мне сейчас важнее разобраться с оценкой точности конечного результата.

Мне к сожалению пока не понятен используемый вами алгоритм: B@R5uk: "Чтобы численно найти частную производную по конкретному аргументу, необходимо немного изменить только один этот аргумент и посмотреть, как изменилась функция." Не могли бы вы дать ссылку с теоретическое обоснование этого подхода и/или объяснить на примере, как я могу перейти от неточности в значении химического потенциала к неточности в итоговой концентрации компонента раствора? Пока я вижу только один путь - это расчитывать минимум при разных исходных значениях $G_0$, изменяющихся в пределах доверительных интервалов. (Боюсь ошибиться, но на мой взгляд, вы предлагаете то же самое.) Чтобы получить статистически значимую выборку таких расчетов должно быть очень много. Пока что во всей мировой литературе по этому вопросу (исключая конечно метод Монте Карло, используя который мне придется расчитывать одно равновесие целый день) мне известна только вот эта книга (стр. 83-121):
https://drive.google.com/file/d/0B08REa ... sp=sharing предлагающая реальное решение. К сожалению, только ее я использовать не могу, потому что мне пока недоступны исходные работы, описывающий алгоритм, и есть подозрение, что в самой книге алгоритм описан неполностью.

[alenov]

Цитата из данной книги, стр. 99:

Цитата:

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

Автор называет этот метод алгоритмом Белова. Интересует адекватное описание этого или подобного алгоритма (лучше, чтобы алгоритм был другой, но с тем же функционалом, и широко известный).

[Deggial]

alenov, оформляйте цитаты тегом quote. Пользуйтесь также кнопками и

[alenov]

Deggial в сообщении #813830 писал(а):

alenov, оформляйте цитаты тегом quote. Пользуйтесь также кнопками и

Deggial, я к сожалению больше не могу редактировать предыдущее сообщение, кончились попытки.

[B@R5uk]

alenov в сообщении #813804 писал(а):

если вы предложите метод, который будет работать быстрее того, который я сейчас использую

Какой метод расчёта вы сейчас используете?

alenov в сообщении #813804 писал(а):

Если константа воды равна 14, то это выражение (ограничение на рН) будет равно 1E-14.

alenov в сообщении #813804 писал(а):

содержание $x_{H^+}$ должно быть около 1E-7 ($pH=-\log_{10} x_{H^+}$

Вы вынуждаете меня разбираться в химических формулах. Водородный показатель -- это минус десятичный логарифм концентрации в моль/литр ионов водорода. В нашей задаче $x_{H^+}$ -- это количество вещества в молях, то есть его надо ещё разделить на объём раствора и только потом логарифмировать. Несколько раз переспрашивая я именно на это и намекал. Итак, в растворе 56,51 молей воды, следовательно её масса 1,018 кг и объём при 25 градусах будет 1,027 литров. ОК, у вас подобрано исходное количество вещества так, чтобы оно численно почти равнялось концентрации. Но так делать нехорошо.

Далее, у меня возникло два вопроса, связанных с формулами, которые я не понимаю, поскольку не химик. Первый вопрос: константа диссоциации воды же дана для чистой воды. А сейчас у нас в системе есть кремний (и в особенности ион $HSiO_3^-$). Не будет ли она (константа) другая в этом случае. И связанный с этим второй вопрос: если рассматривать систему, состоящую исключительно из водорода и кислорода два к одному, то не получится ли из этой формулы как раз водородный показатель 7? Не потому ли:

alenov в сообщении #813804 писал(а):

В литературе по данному вопросу можно найти указание на необходимость условия электронейтральности, про произведение воды (рН) нигде ничего не сказано

что это условие (водородное произведение) для чистой воды получается автоматически если рассматривать состояние системы как минимум функции Гибса-Дюгема?

alenov в сообщении #813804 писал(а):

Еще у вас получилось сильно завышенное содержание кислорода. Реально оно не может быть выше 1Е-60

Ну, у меня получилось, что кислорода в системе $0$ моль с точностью $10^{-16}$ моль. Куда ж ещё точнее?

[ИСН]

alenov в сообщении #813804 писал(а):

Боюсь, что бывают. При учете неидеальности растворов и введении в эту функцию коэффециентов активности

Но Вы-то считаете раствор идеальным, а при этом, как говорится там же, функция выпукла и всё хорошо!

[B@R5uk]

alenov, а откуда вы брали вот эти числа:

alenov в сообщении #812131 писал(а):

Phase_______$G_0_i$($\Delta$), cal/mol
$H_2O,aq$___ -0056690
$H^+$________ 0000000
$H_2,aq$_____ 0004236
$HSiO_3^-$___ -0242801
$O_2,aq$_____ 0003954
$OH^-$______-0037595
$SiO_2,aq$___-0199190

Есть какой-то справочник с ними?

[alenov]

Цитата:

B@R5uk в сообщении #813847 писал(а):

alenov в сообщении #813804 писал(а):

если вы предложите метод, который будет работать быстрее того, который я сейчас использую

Какой метод расчёта вы сейчас используете?

метод сопряженных градиентов Флетчера-Ривса, целевая функция состоит из функции Гиббса-Дюгема, множителей Лагранжа и штрафной функции (в первом сообщении была ссылка)

Цитата:

alenov в сообщении #813804 писал(а):

Если константа воды равна 14, то это выражение (ограничение на рН) будет равно 1E-14.

alenov в сообщении #813804 писал(а):

содержание $x_{H^+}$ должно быть около 1E-7 ($pH=-\log_{10} x_{H^+}$

Вы вынуждаете меня разбираться в химических формулах. Водородный показатель -- это минус десятичный логарифм концентрации в моль/литр ионов водорода. В нашей задаче $x_{H^+}$ -- это количество вещества в молях, то есть его надо ещё разделить на объём раствора и только потом логарифмировать. Несколько раз переспрашивая я именно на это и намекал. Итак, в растворе 56,51 молей воды, следовательно её масса 1,018 кг и объём при 25 градусах будет 1,027 литров. ОК, у вас подобрано исходное количество вещества так, чтобы оно численно почти равнялось концентрации. Но так делать нехорошо.

Не думаю, что вам нужно что-то пересчитывать. Оно и так решается. Хотя действительно, я предложил случайный пример, чтобы получить истинное значение рН, нужно мольное количество ионов водорода помножить на 55.51/53.5102, но это не сильно изменит результат. Спасибо, вы указали мне на мою ошибку, я забыл про это в своем коде!

Цитата:

Далее, у меня возникло два вопроса, связанных с формулами, которые я не понимаю, поскольку не химик. Первый вопрос: константа диссоциации воды же дана для чистой воды. А сейчас у нас в системе есть кремний (и в особенности ион $HSiO_3^-$). Не будет ли она (константа) другая в этом случае.

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

Цитата:

И связанный с этим второй вопрос: если рассматривать систему, состоящую исключительно из водорода и кислорода два к одному, то не получится ли из этой формулы как раз водородный показатель 7? Не потому ли:

alenov в сообщении #813804 писал(а):

В литературе по данному вопросу можно найти указание на необходимость условия электронейтральности, про произведение воды (рН) нигде ничего не сказано

что это условие (водородное произведение) для чистой воды получается автоматически если рассматривать состояние системы как минимум функции Гибса-Дюгема?

Да, возможно именно это и должно работать. Но у меня так не получилось.

Цитата:

alenov в сообщении #813804 писал(а):

Еще у вас получилось сильно завышенное содержание кислорода. Реально оно не может быть выше 1Е-60

Ну, у меня получилось, что кислорода в системе $0$ моль с точностью $10^{-16}$ моль. Куда ж ещё точнее?

1Е-16 в подобных задачах довольно высокая величина. На практике это может приводить к самопроизвольному накоплению кислорода в системе (равновесия сейчас считаются много раз подряд и этот расчетный артефакт обязательно вылезет), где его не может быть и очень серьезным ошибкам. В некоторых задачах важна разница между 1Е-100 и 1Е-200. Но это конечно экзотические случаи, обычно достаточно, чтобы содержание могло быть 1Е-100. Вообще содержание кислорода - важный параметр, одни минералы буферируют его на уровне 1Е-20, другие 1Е-40 и т.д.. В пустой системе (а она фактически пустая) его не может быть больше 1Е-60. Возможность получить значение 1Е-100 одно из требований для подобных алгоритмов.

-- 13.01.2014, 20:50 --

ИСН в сообщении #813851 писал(а):

alenov в сообщении #813804 писал(а):

Боюсь, что бывают. При учете неидеальности растворов и введении в эту функцию коэффециентов активности

Но Вы-то считаете раствор идеальным, а при этом, как говорится там же, функция выпукла и всё хорошо!

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

-- 13.01.2014, 21:07 --

B@R5uk в сообщении #813884 писал(а):

alenov, а откуда вы брали вот эти числа:

alenov в сообщении #812131 писал(а):

Phase_______$G_0_i$($\Delta$), cal/mol
$H_2O,aq$___ -0056690
$H^+$________ 0000000
$H_2,aq$_____ 0004236
$HSiO_3^-$___ -0242801
$O_2,aq$_____ 0003954
$OH^-$______-0037595
$SiO_2,aq$___-0199190

Есть какой-то справочник с ними?

Эти значения рассчитаны в сторонней программе. Есть множество справочников и баз данных для разных целей. Если вам это нужно я могу поделиться тем, что вам нужно, но их очень много. Можно посмотреть здесь: https://drive.google.com/file/d/0B08REa ... sp=sharing Эта статья - классический источник термодинамической информации для подобных задач. А вот одна из последних баз данных такого рода: https://drive.google.com/file/d/0B08REa ... sp=sharing

-- 13.01.2014, 21:10 --

Вы мне лучше про равномерную сетку в гиперкубе чего-нибудь скажите. А то пока мы тут термодинамическим моделированием занимаемся у меня работа не движется совсем:)

[alenov]

Мне посоветовали решить эту задачу следующим очень простым способом. Нахождение неопределенности решений может быть сделано путем решения задачи оптимизации:
$\max f(G') = \sum_i (G^d^e^c_i-G'_i)^2$ при $\left|G^d^e^c_0_i-G'_0_i \right| \leqslant $\Delta$_i$
$G^d^e^c_0_i$ - исходное значение стандартной энергии Гиббса, при котором найдено первоначальное решение, $G'_0_i$ - любое значение этого параметра в эллипсоиде неопределенности (входные данные для задачи оптимизации), $G^d^e^c_i$ - значение свободной энергии при значении $x_i$, найденных в результате решения исходной задачи (в данной задаче это массив констант), $G'_i$ - ищется на каждом шаге при заданных $G'_0_i$. Напомню, что $G'_i = x_i (G'_0_i + \ln (x_i/\sum_i x_i))$.
Я так понимаю, что изначально предложенное решение B@R5uk'а по сути является частным случаем этой задачи, если предположить, что многомерному эллипсоиду $G'_0_i$ соответствует аналогичный эллипсоид $x_i$, что к сожалению почти всегда не соответствует действительности.
У меня вопрос к уважаемым участникам форума, как эту задачу оптимизации решить за наименьшее количество шагов (вычисления $x_i$ на каждом шаге требуют много времени)? Я могу применить максимизацию с штрафными функциями, но мне так кажется, что это далеко не самый оптимальный подход, т.к. целевая функция выглядит очень простой. Точность найденного решения в данном случае значения практически не имеет, лишь бы можно было математически обоснованно сказать, что найденные $G'_i$ действительно соответствуют максимуму $f(G')$.

[alenov]

На данный момент я остановился на следующем методе. Если точка максимума не попадает в область возможных значений, то она очевидно лежит в вершине этого ящика. Я нахожу максимальную длину вектора за 2 шага. На первом считаю значение целевой функции при всех максимальных $G'_0_i$, рассчитываю градиент в этой точке, если по каким-то направлениям градиент отрицательный, то для этих направлений беру минимальные $G'_0_i$, считаю значение целевой функции в этой точке. Если знаки градиента не изменились, значит это и есть максимум. Мой пример вроде бы отвечает этой схеме. Однако возможна ситуация, когда максимум внутри ящика. Что делать в этой ситуации, которая возможна, пока не решил. По-видимому, нужно убедиться, перебрав все вершины, что она внутри, а затем искать максимум методом сопряженных градиентов с одномерным поиском методом золотого сечения.