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

[alenov]

Передо мной стоит задача оценки точности/достоверности найденного решения задачи по минимизации функции. Алгоритм находит значения, которые соответствуют минимальному значению функции при заданных ограничениях, затем нужно найти что-то вроде сигмы или доверительного интервала этих значений (для каждого в отдельности). Я пока очень слабо представляю с чего начать. Краем уха слышал, что подобную задачу можно решить с помощью метода Монте Карло. Также для этого можно использовать рассчитанные значения множителей Лагранжа. Но кроме краткого упоминания в Википедии и в одном из учебников я пока не нашел ни одного примера как это делается. Физический смысл у этой задачи следующий. Например, рассчитывается равновесный состав газов в реакторе. В реальном реакторе концентрация газов всегда определяется с некоторой точностью, обусловленной методикой определения и отбора пробы. Вот мне как раз и нужно расчитать этот модельный аналог точности, обусловленной как точностью исходных параметров (энергия Гиббса каждого компонента, их активность и ограничения), так и точностью расчета минимума функции. В каком направлении можно двигаться?

Я занимался поиском минимума функции в этой теме:
topic76893.html С вашей помощью эта задача решена. Хочу еще раз поблагодарить Мат Ламера за ряд полезных рекоммендаций.

[B@R5uk]

Если про погрешности прямого измерения вы знаете, то уж точно должны знать про погрешности косвенного измерения. Если не знаете, то первое, что мне приходит на ум это книжка Митин, Русаков "Анализ и обработка экспериментальных данных" Более новая версия более простая, только для ознакомления. У Митина есть ещё старая книжка, в которой рассказано как считать погрешность косвенного измерения, найденного нелинейным МНК численно.

Конкретно в вашем случае я не знаю, что надо делать, так как не вижу конкретной формулы.

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

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

$x=\frac{1}{2}(x'+\sqrt{x'^2+4})=\frac{2}{-x'+\sqrt{x'^2+4}}$

Эта функция делает ту же работу, но ведёт себя более спокойно.

[alenov]

Спасибо, B@R5uk! Целая куча полезной информации. Начну с книги Митина, Русакова.

А вот предложенное преобразование я, к сожалению, не очень понимаю. Что такое $x'$ и $x$ - исходная и синтетическая переменные? Как же тогда может выполняться равенство? Или первый случай нужно применять, если $x' \geq 0$, а второй случай если $x' < 0$, наверное так. Попробую.

[B@R5uk]

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

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

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

-- 08.01.2014, 18:20 --

Вопрос. Сколько у вас условий вида

$\sum\limits_{i=1}^N(a_{ij}x_i)=b_j$

То есть сколько значений пробегает параметр $j$?

Ещё вопрос. При практических расчётах часто ли бывает так, что переменная $x_i$ будет оптимальная, если будет меньше нуля, но приходится воспользоваться условием $x_i>0$?

[alenov]

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

$\min G(x)=\sum(x_i(G_0_i+RT\ln(x_i/\sum(x_i))))$, $G_0_i = \operatorname{const}, R = 1.985877534 cal/°K/mol, T = 298.15°K$

В этом выражении неточность измерения содержат только $G_0_i$. R и T можно считать абсолютно точными. Задачей же является расчитать неточность для $x_i$.

i может быть порядка 100, j - порядка 20. Заменой $ x = \exp(x')$ я добился, что x не может быть отрицательным. Насколько я понимаю, предложенная вами операция делает то же самое. Если бы этого не было, то конечно, на первой же иттерации какой-то из $x_i$ становится меньше нуля. Оптимальный $x_i$ никак не может быть меньше нуля, потому что это значение концентрации.

Ниже пример найденного минимума, чтобы задача была понятнее.

Phase_______$G_0_i$($\Delta$), cal/mol____$x_i$($\Delta$), mol

$H_2O,aq$___-0056690(5)_________53.509999812418(?)
$H^+$________0000000(5)_________1.87716779053972E-7(?)
$H_2,aq$_____0004236(5)_________1.99999999993131(?)
$HSiO_3^-$____-0242801(5)________1.34312071062245E-7(?)
$O_2,aq$_____0003954(5)_________2.09002173538838E-88(?)
$OH^-$______-0037595(5)________5.32717429437457E-8(?)
$SiO_2,aq$___-0199190(5)________9.98656867668485E-5(?)

$G(x) = -3030128.76598301$
В данном случае $i = 7$, а $j = 5$ (O, H, Si, z(эллектронейтральность) и pH). Истинное значение доверительного интервала ($\Delta$) для $G_0_i$ надо еще уточнить, пусть пока будет равет 5.

[B@R5uk]

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

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

Нет. В книжке МНК приведён для простейшего случая, а именно прямой зависимости. Вообще говоря, МНК можно спокойно применять для любой сколь угодно нелинейной зависимости. К вашей проблеме МНК не применим под другой причине.

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

У вас совсем другая ситуация. У вас есть набор измеренных величин и вы по ним рассчитываете другой набор величин. Исходные величины имеют погрешность. Следовательно, рассчитанные тоже имеют погрешность.

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

Ниже пример найденного минимума, чтобы задача была понятнее.

Да, ситуация для меня немного проясняется. По табличке я вижу, что вы ищете концентрации семи веществ. А условия

$\sum\limits_{i=1}^N(a_{ij}x_i)=b_j$

на них какие-то были наложены?

[ИСН]

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

А условия $\sum\limits_{i=1}^N(a_{ij}x_i)=b_j$ на них какие-то были наложены?

Валовый состав и электронейтральность, очевидно.

[B@R5uk]

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

Валовый состав и электронейтральность, очевидно.

Мне нужны конкретные цифры, чтобы я во-первых сам до конца разобрался откуда берётся решение, во-вторых проверил своим расчётом результат (на всякий случай), в-третьих на конкретном примере объяснил, как считать погрешность.

-- 09.01.2014, 22:57 --

В общем случае, когда есть функция $F(x_1,x_2,...,x_N,c_1,c_2,...,c_M)$, минимум которой надо найти и которая зависит от искомых переменных $x_k$ и от заданных параметров $c_l$. Если набор переменных $x_k$ доставляет минимум функции $F$, то они удовлетворяют системе из $N$ уравнений:

$\frac{\partial F}{\partial x_n}(x_1,x_2,...,x_N,c_1,c_2,...,c_M)=0$

Пусть решение этой системы уравнений -- функции вида (в общем случае не записываемые явно)

$x_k=X_k(c_1,c_2,...,c_M)$

Тогда погрешности величин $x_k$ будут выражаться через погрешности величин $c_m$ и производные

$\frac{\partial X_k}{\partial c_m}(c_1,c_2,...,c_M)$

которые находятся по теореме о неявной функции из системы уравнений

$\frac{\partial F}{\partial x_n}(x_1,x_2,...,x_N,c_1,c_2,...,c_M)=0$

Но это простой случай, где ищется обычный минимум, у вас в задаче минимум условный, поэтому необходимые для расчёта погрешностей производные функций $X_k$ будут рассчитываться по-другому. Пока пытаюсь сообразить как.

[alenov]

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

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

Валовый состав и электронейтральность, очевидно.

Мне нужны конкретные цифры, чтобы я во-первых сам до конца разобрался откуда берётся решение, во-вторых проверил своим расчётом результат (на всякий случай), в-третьих на конкретном примере объяснил, как считать погрешность.

O - 53.5102 (cостав системы в моль по кислороду, водороду и кремнию)
H - 111.02
Si - 0.0001
z - 0 (отрицательные и положительные заряды должны компенсировать друг друга)
pH ($-\lg(H^+)-\lg(OH^-)$ ) - 14 (ну это тоже нужно, долго рассказывать зачем, для понимания задачи не принципиально, можно решить и без этого условия, но с ним получается точнее)

-- 09.01.2014, 23:30 --

Для кислорода это ограничение будет такое:
$b_O = 1 \cdot x_H_2_O + 0 \cdot x_H_+ + 0 \cdot x_H_2 + 3 \cdot x_H_S_i_O_3_- + 2 \cdot x_O_2 + 1 \cdot x_O_H_- + 2 \cdot x_S_i_O_2$ = 53.5102 mol

-- 09.01.2014, 23:36 --

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

$\frac{\partial X_k}{\partial c_m}(c_1,c_2,...,c_M)$

Похоже на множители Лагранжа. Разве нет?

[B@R5uk]

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

Похоже на множители Лагранжа. Разве нет?

Нет, это совсем другое.

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

Заменой $ x = \exp(x')$ я добился, что x не может быть отрицательным. Насколько я понимаю, предложенная вами операция делает то же самое.

Да, совершенно верно.

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

pH ($-\lg(H^+)-\lg(OH^-)$ ) - 14

К сожалению, не понял, как эту связь выразить в формулу.

С трудом получил те же самые числа (здесь $y$ -- это переменные, по которым велась оптимизация без учёта pH):

>> G
G =
   -3.030128769625446e+006
>> x
x =
  53.509993249495970
   0.000006750504023
   2.000000000000007
   0.000002412632721
   0.000000000000004
   0.000004337871302
   0.000097587367279
>> y
y =
   0.999999873845935
   0.000000126154065
   0.024126327209463

С трудом -- потому что решение "в лоб", которое легко изменяется под любую другую такую же задачу и легко масштабируется на любое другое число переменных/связей, не заработало. Пришлось придумать другое решение, красивое, но заточенное под эту конкретную задачу. Масштабировать под другую задачу его тоже можно, то только в ручную, автоматически это сделать будет затруднительно. Если вас интересует могу выложить исходники.

Впрочем, для расчёта погрешностей не важно каким способом получено решение. Формулы одни и те же.

Кстати, простейший вариант расчёта погрешностей вычисляемых величин -- это изменять по одной каждую из величин на величину погрешности в ту и другую сторону и пересчитывать искомые иксы. Фактически это будет численным вычислением производной.

-- 10.01.2014, 20:52 --

Объясните, пожалуйста, откуда берутся эти числа:

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

O - 53.5102
H - 111.02
Si - 0.0001

Это измеряемые величины (тогда у них тоже есть погрешность) или же это те величины, которые надо получить? Кстати, это количество вещества в молях, или же концентрация в молях на литр?

[alenov]

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

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

Да, предложенный вами метод мне тоже приходил в голову, но там нужно попробовать комбинацию всех крайних (максимальных - минимальных) и средних значений со всеми, т.е. $7 \cdot 7 \cdot 7$. Это уже довольно много расчетов получается. А эта система существенно упрощена, в нормальной задаче будет порядка сотни компонентов, это уже миллион расчетов. Я нашел указание, что подобные задачи можно решать с помощью "гиперкуба", состоящего из равноудаленных точек. С помощью этого подхода количество проб в задаче с размерностью 100 входных элементов с погрешностью можно снизить до 100-1000. Пока нормального описания этого метода не нашел.

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

[ИСН]

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

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

Исключено. В химии не бывает лишних уравнений, т.е. переопределённых систем. И недоопределённых тоже.
Кстати, а Вы учитываете возможность разделения фаз? Типа что-то взлетело, или что-то выпало в осадок?

[alenov]

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

[B@R5uk]

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

Собираюсь учитывать, пока рассматриваю простейший случай. Нужно отработать алгоритм.

Ну, это не верный подход. Та задача, которую вы сейчас поставили очень хорошая: целевая функция выпукла (если не ошибаюсь) на множестве определения и задана так же на выпуклом множестве. В таком хорошем случае минимум существует и единственен. То есть если программа находит какой-то локальный минимум (с учётом ограничений), то это будет как раз тот самый, который нужен. Если же добавить нелинейное ограничение $(x_{OH^-})\cdot(x_{H^+})=\operatorname{const}$ (кстати, чему равна эта константа, если ${x}_{k}$ -- количество вещества в молях?), то даже выпуклая функция может иметь несколько локальных минимумов. В этом случае программа должна найти их все и выбрать наименьший, а это весьма сложная задача.

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

Также интересует, за сколько времени рассчитывается полученное вами решение.

Время вычислений с точностью до $10^{-16}$ (максимальная возможная точность) доли секунды и это учитывая то, что программа на матлабе (в том числе и функция многомерной оптимизации) -- это скрипт, а не код. Для достижения указанной точности при указанных вами значениях констант и связей целевая функция вычислялась 1060 раз.

Использовать вашу подстановку $x=\exp(x')$ или мою с иррациональностью нельзя, так как в случае, когда минимум функции достигается, когда один из иксов равняется нулю, соответствующий штрихованный икс должен равняться минус бесконечности, то есть после такой подстановки целевая функция получается не имеет экстремума (он лежит на бесконечности). Этим объясняется то, что ваша программа очень долго считает: одна из координат долго уходит на бесконечность (так как делает это конечными шагами).

Чтобы избежать эту проблему, но, тем не менее, удовлетворить условиям положительности иксов, целевую функцию я определил следующим образом:
$G=\sum\limits_{k=1}^{N}{{{x}_{k}}\left( {{G}_{0k}}+RT\ln \frac{{{x}_{k}}}{\sum\limits_{l=1}^{N}{{{x}_{l}}}} \right)}$, если все ${x}_{k}$ положительны
$G=\infty$, если хотя бы один ${x}_{k}$ отрицателен.

Для поиска минимума целевой функции я использовал встроенную в МАТЛАБ процедуру fminsearch, которая использует алгоритм Нелдера — Мида многомерной оптимизации. Этот алгоритм очень удобен, но имеет свои практические недостатки. Например, все аргументы функции должны иметь один порядок. Для обхода этого недостатка, для уменьшения количества переменных, по которым будет вестись оптимизация, а так же для удовлетворения связей, наложенных на аргументы, я использовал следующую подстановку:

$\[\left( \begin{matrix} {{x}_{1}} \\ {{x}_{2}} \\ {{x}_{3}} \\ {{x}_{4}} \\ {{x}_{5}} \\ {{x}_{6}} \\ {{x}_{7}} \\ \end{matrix} \right)=\left( \begin{matrix} 0 \\ 0 \\ 55,51 \\ 0 \\ 26,755 \\ 0 \\ 0,0001 \\ \end{matrix} \right)+\left( \begin{matrix} 53,51 & 0 & 0 \\ 0 & 53,51 & 0 \\ -53,51 & -53,51 & 0 \\ 0 & 0 & 0,0001 \\ -\text{26}\text{,755} & -\text{26}\text{,755} & 0 \\ 0 & 53,51 & -0,0001 \\ 0 & 0 & -0,0001 \\ \end{matrix} \right)\left( \begin{matrix} {{y}_{1}} \\ {{y}_{2}} \\ {{y}_{3}} \\ \end{matrix} \right)\] $

Чтобы получить эту подстановку я исследовал связи и ограничения (благо, что все связи линейные). Нашёл, в каких диапазонах может изменяться каждый из иксов, выбрал зависимые и независимые переменные, выразил зависимые через независимые и независимые отмасштабировал с учётом их диапазона так, чтобы они изменялись в пределах от нуля до единицы. Таким образом, я уменьшил размерность задачи, и построил "хорошую" для процедуры оптимизации функцию.

И так, код программы:

(Оффтоп)

Файл "prov3.m"

clc
clear
format compact
RT=1.985877534*298.15;
G0=[
   -0056690     %   H2O
    0000000     %   H+
    0004236     %   H2
   -0242801     %   HSiO3
    0003954     %   O2
   -0037595     %   OH-
   -0199190];   %   SiO2
G0dRT=G0/RT;
c1=53.51;
c2=26.755;
c3=0.0001;
A=[
    c1  0   0
    0   c1  0
   -c1 -c1  0
    0   0   c3
   -c2 -c2  0
    0   c1 -c3
    0   0  -c3];
B=[
    0
    0
    55.51
    0
    c2
    0
    c3];
fmso=optimset('fminsearch');
fmso=optimset(fmso,'Display','iter','TolFun',1e-20,'TolX',1e-20,'MaxIter',200000,'MaxFunEvals',400000);
y=fminsearch(@(y)prov3_fnc(B+A*y,G0dRT),0.1*ones(3,1),fmso);
x=B+A*y;
G=RT*prov3_fnc(x,G0dRT);
disp('x=')
disp(x)
disp('G=')
disp(G)

Файл "prov3_fnc.m"

function G=prov3_fnc(x,G0dRT)
    if sum(x<0)>0
        G=Inf;
    else
        G=sum(x.*(G0dRT+log(x/sum(x))));
    end
end

-- 13.01.2014, 14:55 --

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

Я не знаю, насколько необходимо это условие.

Это условие влияет на концентрации и, следовательно, на количества вещества ионов $H^+$ и $OH^-$. Эти количества вещества связаны через электронейтральность с количеством вещества ионов $HSiO_3^-$, следовательно, если в растворе не будет достаточно много других заряженных ионов, то введение ограничения, связанного с ионным произведением воды, может изменить концентрацию $HSiO_3^-$ на порядок и более.

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

но там нужно попробовать комбинацию всех крайних (максимальных - минимальных) и средних значений со всеми, т.е. $7 \cdot 7 \cdot 7$

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

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

П.С. И всё-таки, откуда взялась функция
$G=\sum\limits_{k=1}^{N}{{{x}_{k}}\left( {{G}_{0k}}+RT\ln \frac{{{x}_{k}}}{\sum\limits_{l=1}^{N}{{{x}_{l}}}} \right)}$
и почему в ней иксы -- это количество вещества, а не концентрации? Это же раствор.

[ИСН]

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