Потенциометрия — справка ElectroChemLab

1.Назначение

Из многих методов физико-химического исследования свойств растворов электролитов потенциометрический метод является наиболее точным, термодинамически строгим и характеризуется однозначностью вследствие селективности. Раздел программного комплекса «ElectroChemLab» — «Потенциометрия» — предназначен для обработки экспериментальных данных потенциометрических измерений, полученных методом стандартных серий: измерений ЭДС гальванических элементов с переносом и без переноса в растворах со специальным образом подобранным составом.

По измеренным значениям ЭДС гальванического элемента производится итерационное уточнение (оптимизация) оценок выбранных параметров химической модели раствора — стандартной ЭДС гальванического элемента $E^{\circ}$, логарифмов термодинамических констант равновесий $\lg K_i\ (i = 0 \div s)$, параметров расширенного уравнения Дебая — Хюккеля для коэффициентов активности ионов ($\mathring{a}$, $b$). Практика показала, что уравнение такого вида описывает коэффициенты активности ионов в широком диапазоне ионных сил с помощью одного эмпирического параметра в растворах средних концентраций.

Расчёт равновесного состава для каждого раствора выполняется по универсальному алгоритму, представленному в разделе «Равновесия» (матричный метод описания равновесий, основанный на формализме Бринкли; итерационный процесс нахождения корней системы нелинейных уравнений по методу Ньютона). Решением прямой задачи химических равновесий являются логарифмы активности частиц базиса $\ln a_r\ (r = 1 \div q)$; исходными параметрами являются значения термодинамических констант равновесия и параметры расширенного уравнения Дебая — Хюккеля. В процессе уточнения логарифмов активности частиц базиса происходит и итеративное уточнение ионной силы раствора.

См. также: справку раздела «Расчёт равновесного состава» — там подробно описаны матричный формализм Бринкли, расширенное уравнение Дебая — Хюккеля и алгоритм решения прямой задачи, которая здесь используется как «чёрный ящик».

2.Идеология поиска решения

Процесс поиска адекватной модели химических равновесий и её количественного описания сводится к минимизации остаточной суммы квадратов отклонений значений электродвижущей силы (ЭДС), рассчитанных по математической модели эксперимента, и прямых экспериментальных данных по методу наименьших квадратов. Подбираются такие значения вектора параметров $\vec{\theta} = (E^{\circ},\ \lg K_i\ (i = 0 \div s),\ \mathring{a},\ b)$, при которых рассчитанные модельные значения ЭДС наилучшим образом (в смысле взвешенной суммы квадратов) описывают измеренные величины на всём интервале ионных сил.

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

Двухуровневая схема

Внешний уровень — оптимизатор Варьирует вектор выбранных параметров (например, $E^{\circ}$ и несколько $\lg K$), решая обратную задачу химических равновесий — математическую задачу определения параметров системы на основе измеренных в эксперименте равновесных концентраций, связанных функционально с массивом экспериментальных данных.

Внутренний уровень — моделирование Для каждого пробного $\vec{\theta}$ решается прямая задача химического равновесия во всех растворах серии; из равновесных концентраций и коэффициентов активности потенциалопределяющих ионов вычисляется активность, затем — расчётная ЭДС по Нернсту, и наконец — целевая функция.

Разнесение двух источников нелинейности (само химическое равновесие и зависимость искомых величин от параметров) позволяет переиспользовать надёжный равновесный решатель «как есть»: оптимизатор обращается к нему как к «чёрному ящику», получая равновесный состав.

Кэширование. Если варьируются только параметры, не влияющие на равновесие (то есть лишь $E^{\circ}$), равновесие рассчитывается один раз и кэшируется — это существенно ускоряет подбор. При изменении $\lg K$, $\mathring{a}$ или параметра $b$ равновесие пересчитывается.

3.Уравнение Нернста и типы цепей

ЭДС гальванического элемента (в мВ) вычисляется по термодинамически строгому уравнению Нернста:

$$ E_{\text{расч}} \;=\; E^{\circ} \;-\; \frac{RT}{nF}\,\ln \prod_{j}^{p} a_j^{\,\nu_j} \;=\; E^{\circ} \;-\; \frac{\theta}{n} \sum_{j}^{p} \nu_j \lg a_j\,, \qquad a_j = c_j \gamma_j\,, $$

где $n$ — число электронов в электродной реакции, $\nu_j$ — показатель степени при активности $j$-го потенциалопределяющего иона, $a_j$ — активность иона в растворе (произведение равновесной концентрации $c_j$ на коэффициент активности $\gamma_j$, взятые из расчёта равновесия), $\theta$ — коэффициент в уравнении Нернста:

$$ \theta \;=\; \ln(10)\,\frac{RT}{F}\,1000\,; $$

Множитель 1000 переводит вольты в милливольты; при 298,15 K величина $\theta \approx 59{,}16$ мВ на единицу $\lg a$.

Знак логарифмического слагаемого. Знак «−» соответствует электроду, обратимому к катионам (например, водородный электрод), как отрицательному. При использовании ионоселективных электродов возможна замена знака на противоположный — для этого в окне ввода предусмотрен переключатель знака.

Типы цепей

Цепь без переноса ЭДС зависит от активностей двух ионов (катиона и аниона), то есть от среднеионной активности электролита. Поэтому задаются два потенциалопределяющих иона из совокупности частиц в растворе со своими показателями степеней $\nu$ в уравнении Нернста.

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

ЭДС гальванического элемента и $E^{\circ}$ задаются и обрабатываются в милливольтах, так как при этом уровень погрешности сопоставим по величине с погрешностями остальных оптимизируемых параметров $\vec{\theta}^{\,T} = (\lg K_i\ (i = 0 \div s),\ \mathring{a},\ b)$, что удобно для совместной оптимизации.

4.Целевая функция и веса

Рис. 1. Этап 5 «ЭДС»: тип цепи, число электронов, потенциалопределяющие ионы и таблица измерений. Концентрации в таблице только для чтения — они приходят с этапов ввода равновесия, а вводится ЭДС и, при нескольких повторностях, веса. Пример — серия ацетатных буферных растворов со стеклянным электродом: по ней определяются одновременно E° ячейки и константа кислотности.

Целевая функция $F(\vec{\theta})$ — это математическое выражение, которое подлежит минимизации посредством оптимизации параметров $\vec{\theta}$ для получения точного решения задачи. Целевая функция определяет желаемый результат в смысле метода наименьших квадратов — истинное решение, при котором достигается максимальное совпадение расчётных и экспериментальных значений ЭДС. Это достигается поиском минимума остаточной суммы квадратов отклонений (RSS), построенной на невязках $\varepsilon_k = E_{\text{расч},k} - E_{\text{эксп},k}$.

При проведении эксперимента измеренные значения ЭДС всегда содержат погрешность, которая может быть различной. При точных исследованиях каждый раствор исследуется несколько раз ($m$) при совокупном определении параметров, и используют наиболее вероятные значения, которыми являются средние арифметические $\overline{E}_k$. Для учёта «вклада» конкретного измерения в итоговые результаты используют веса — множители, обратно пропорциональные погрешностям, что даёт меньший вклад в конечный результат менее точным измерениям. Веса, как правило, нормируются. Если же измерения гомогенны (выполнены с одинаковой погрешностью), то весам назначается значение, равное 1.

В приложении реализован наиболее общий случай поиска адекватного решения, основанный на минимизации взвешенной остаточной суммы квадратов отклонений:

$$ F(\vec{\theta}) \;=\; \sum_{k}^{N} w_k\,\varepsilon_k^{2} \;=\; \sum_{k}^{N} w_k \left( E_{\text{расч},k} - E_{\text{эксп},k} \right)^{2}. $$

Суммирование взвешенных невязок осуществляется по всем $N$ растворам серии; $\overline{E}_{\text{эксп},k}$ — среднее по повторным измерениям в $k$-ом растворе.

Весовая матрица

По умолчанию все веса $w_k = 1$. Если в растворе выполнено более одного измерения, вес вычисляется как обратная выборочная дисперсия:

$$ w_k \;=\; \frac{1}{s_k^{2}}\,; \qquad\qquad s_k^{2} \;=\; \frac{1}{m_k - 1} \sum_{j}^{m_k} \left( E_{kj} - \overline{E}_k \right)^{2}. $$

Явно заданные пользователем веса имеют приоритет над расчётными. Веса вводятся на этапе ввода концентраций и ЭДС.

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

5.Вектор параметров

Оптимизатор работает с единым вектором $\vec{\theta}$, который собирается из разнородных по физико-химическому смыслу величин:

Параметр	Физико-химический смысл
$E^{\circ}$	стандартная ЭДС элемента (может быть известна и зафиксирована, либо уточняться как параметр)
$\lg K$	логарифмы констант равновесий; выбираются не все, а только нужные — по стехиометрической матрице, в соответствии с равновесием образования конкретной частицы
$\mathring{a}$	параметр максимального сближения ионов
$b$	эмпирический параметр уравнения для коэффициентов активности

Для каждого параметра — предположительного кандидата на уточнение — задаётся, независимо от вовлечения в процесс оптимизации: начальное приближение и допустимые границы (минимальное и максимальное значения). Границы в процессе оптимизации учитываются мягким штрафом (безградиентные методы формально не ограничены).

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

Набор параметров зависит от выбранной модели активности

Параметры $\mathring{a}$ и $b$ в списке кандидатов появляются только тогда, когда их читает выбранная на шаге 1 модель коэффициентов активности — то есть при расширенном уравнении Дебая–Хюккеля. У предельного закона и уравнения Дэвиса подгоняемых параметров нет вовсе, у MSA единственный параметр — радиус иона.

Это не косметическое ограничение интерфейса. Подгонка по параметру, которого модель не читает, двигала бы величину, ни на что не влияющую: соответствующая строка и столбец матрицы Гессе оказались бы нулевыми, гессиан — вырожденным, а доверительные интервалы всех параметров — недостоверными. Пользователь получил бы предупреждение о вырожденности вместо объяснения, поэтому такие строки просто не предлагаются.

Методология коэффициентов активности изложена в справке «Расчёт равновесного состава», раздел 5. Там же — лестница из четырёх моделей с уравнениями и границами применимости, разбор термодинамической согласованности, порядок выбора модели и то, что выбор модели меняет в потенциометрии (раздел 5.7). Эта страница описывает только поиск решения.

6.Методы оптимизации

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

По желанию пользователя при уже найденном решении (близком к точке минимума функции) всегда можно уточнить координаты точки решения подключением на завершающем этапе метода Ньютона второго порядка. Этот метод использует матрицу вторых производных — матрицу Гессе $\mathbf{H}$, которая несёт информацию о кривизне поверхности функции в пространстве параметров. Эта особенность заключения расчёта отмечается «флажком».

В программном комплексе представлены следующие методы оптимизации.

6.1. Метод Нелдера — Мида (Nelder — Mead)безградиентный

Также известен как метод деформируемого многогранника и симплекс-метод^[1] — метод безусловной оптимизации функции многих переменных, не использующий производной (градиента) функции. Суть метода заключается в последовательном перемещении и деформировании симплекса вокруг точки экстремума. Метод устойчив к шуму и негладкости целевой функции. Эффективен при небольшом числе параметров (примерно до 5–10); с ростом размерности сходимость замедляется.

Операции: отражение, расширение, сжатие (внешнее/внутреннее), редукция симплекса. Использован двойной критерий остановки оптимизации — по разбросу значений функции на симплексе и по размеру симплекса.

6.2. Метод покоординатного спуска Хука — Дживса (Hooke — Jeeves)безградиентный

Также известен как метод конфигураций^[2]. Как и алгоритм Нелдера — Мида, служит для поиска безусловного локального экстремума функции и относится к прямым методам, то есть опирается непосредственно на значения функции. Он работает по принципу «шаг за шагом». Алгоритм делится на две фазы, чередующиеся вокруг базовой точки:

Исследующий поиск — по очереди меняется каждый параметр, чтобы выяснить, куда ведёт наклон: покоординатно пробуется шаг $\pm\delta_i$, принимается любое улучшение (нащупывается спуск).
Образцовый (pattern) ход — «прыжок» вперёд в том направлении, где нашли спуск: экстраполяция в направлении суммарного сдвига $\mathbf{x}_{new} = \mathbf{x} + (\mathbf{x} - \mathbf{x}_{base})$, «разгоняющая» спуск вдоль найденного оврага (компенсирует покоординатность).

Если ни исследование, ни образец не улучшают функцию — шаг дробится ($\delta \leftarrow \rho\,\delta$, $\rho = 0{,}5$), и поиск продолжается на более мелкой сетке, пока не будет найдена точная точка. Начальный шаг $\delta_i = |x_i^{\circ}|\,\rho$ сам подстраивается под масштаб каждого параметра — это важно для смеси разномасштабных величин (например, $E^{\circ} \approx 250$ мВ и $\lg K \approx 4$).

Метод нагляден и полезен как быстрый и устойчивый предварительный проход при небольшом числе параметров и при слабой связанности параметров; на сильно «овражных» функциях медленнее методов с сопряжёнными направлениями. Границы учитываются тем же штрафом; критерий останова — по относительному размеру шага.

6.3. Метод Левенберга — Марквардта (Levenberg — Marquardt)градиентный

Метод оптимизации для решения задач о наименьших квадратах^[3]. Представляет собой комбинацию метода Ньютона и метода градиентного спуска — демпфированный МНК (диагональное приближение гессиана). Даёт быструю сходимость на гладкой целевой функции.

Так как целевая функция возвращает RSS, градиент и диагональ гессиана вычисляются конечными разностями, а шаг

$$ \delta_i \;=\; -\,\frac{g_i}{H_{ii} + \lambda} $$

демпфируется параметром $\lambda$: он уменьшается при успешном шаге (ньютоновское поведение вблизи минимума) и растёт при неудаче (скатывание к градиентному спуску). Устойчив при плохой обусловленности (высокая ионная сила, корреляция параметров). Границы учитываются тем же штрафом, что и для Нелдера — Мида. Критерии останова — по норме градиента и по размеру шага.

6.4. Метод Гаусса — Ньютона (Gauss — Newton)градиентный

Используется для решения задач нелинейным методом наименьших квадратов^[4]. Является модификацией метода Ньютона, но, в отличие от него, может быть использован только для минимизации суммы квадратов. Его преимущество в том, что не требуется вычисления вторых производных. Теоретически самый быстрый метод на гладких задачах МНК — вблизи минимума сходимость близка к квадратичной.

В отличие от диагонального Левенберга — Марквардта (который работает прямо по скалярной RSS), метод использует явный вектор невязок

$$ r_k \;=\; \sqrt{w_k}\,\left( E_{\text{расч},k} - E_{\text{эксп},k} \right) $$

и строит полный якобиан $\mathbf{J} = \partial \mathbf{r} / \partial \vec{\theta}$ (центральными разностями), а аппроксимация гессиана $\mathbf{J}^{T}\mathbf{J}$ учитывает связь параметров (недиагональные элементы). Направление $\delta$ решает нормальные уравнения:

$$ (\mathbf{J}^{T}\mathbf{J})\,\delta \;=\; -\,\mathbf{J}^{T}\mathbf{r}. $$

Оптимальный выбор шага («демпфированный Гаусс — Ньютон»): вдоль $\delta$ выполняется точный одномерный поиск методом Брента по RSS — множитель $\alpha$ подбирается оптимально. Вблизи решения $\alpha \to 1$ (полный ньютоновский шаг, быстрая сходимость), вдали шаг демпфируется, что предотвращает расходимость «чистого» Гаусса — Ньютона на нелинейных участках — без введения параметра Марквардта $\lambda$. При вырождении $\mathbf{J}^{T}\mathbf{J}$ шаг откатывается к наискорейшему спуску.

Границы учитываются тем же штрафом; критерии останова — по норме $\|\mathbf{J}^{T}\mathbf{r}\|$ и по размеру шага. Наиболее эффективен на гладких задачах с малым — средним числом параметров.

6.5. Метод BFGS (Broyden — Fletcher — Goldfarb — Shanno)квазиньютоновский

Итерационный метод численной оптимизации, предназначенный для нахождения локального экстремума нелинейного функционала без ограничений^[5]. BFGS — один из наиболее широко применяемых квазиньютоновских методов, в которых не вычисляется напрямую гессиан функции: гессиан оценивается приближённо, исходя из сделанных до этого шагов.

Метод хранит и обновляет явную аппроксимацию обратного гессиана $\mathbf{H}$ (ранг-2 формула BFGS), задаёт направление поиска $\mathbf{p} = -\mathbf{H}\mathbf{g}$ и выполняет вдоль него точный одномерный поиск методом Брента. Даёт быструю сверхлинейную сходимость на гладкой целевой функции и учитывает связь параметров (в отличие от диагонального приближения Левенберга — Марквардта), но требует памяти $O(n^2)$ на матрицу $\mathbf{H}$. При нарушении условия кривизны ($\mathbf{y}^{T}\mathbf{s} \le 0$) матрица $\mathbf{H}$ перезапускается в единичную.

Границы учитываются тем же штрафом; критерии останова — по норме градиента и по размеру шага. Данный метод находит минимум любой дважды непрерывно дифференцируемой выпуклой функции. Несмотря на эти теоретические ограничения, как показывает опыт, BFGS хорошо справляется и с невыпуклыми функциями.

6.6. Метод Пауэлла (Powell) — сопряжённые направлениябезградиентный

Не требует производных^[6]. За внешнюю итерацию последовательно минимизирует функцию по $n$ направлениям точным одномерным поиском Брента. Затем вдоль нового составного направления

$$ \mathbf{d} \;=\; \mathbf{x}_e - \mathbf{x}_s $$

(суммарный сдвиг от начала к концу итерации) выполняется дополнительный линейный поиск, и это направление заменяет в наборе направление наибольшего убывания — так постепенно строятся приближённо сопряжённые направления, ускоряющие спуск в «овражных» задачах. Обновление набора защищено тестом Флетчера (классическая форма), предотвращающим вырождение (линейную зависимость) направлений.

Эффективен при небольшом числе параметров (до 15). Границы учитываются тем же штрафом; критерии останова — по относительному изменению функции и по размеру шага.

6.7. Метод сопряжённых градиентов (Polak — Ribière+)градиентный

Высокоэффективный алгоритм сопряжённых градиентов для поиска локального минимума многомерных функций^[7], превосходящий классический метод Флетчера — Ривса благодаря лучшей сходимости на неквадратичных функциях и меньшей чувствительности к ошибкам округления. Градиентный метод, реализованный без хранения матрицы в памяти. Направление

$$ \mathbf{d}_{k+1} \;=\; -\mathbf{g}_{k+1} + \beta\,\mathbf{d}_k \,, \qquad \beta \;=\; \max\left( 0,\; \frac{\mathbf{g}_{k+1}(\mathbf{g}_{k+1} - \mathbf{g}_k)}{\|\mathbf{g}_k\|^{2}} \right) $$

с коэффициентом Полака — Рибьера (вариант PR+ обрезает отрицательные $\beta$, гарантируя направление спуска); вдоль направления — точный линейный поиск Брента, градиент — конечными разностями.

Главное преимущество реализации — не хранит матрицу: $O(n)$ памяти против $O(n^2)$ у BFGS, поэтому выгоден при большом числе параметров. Направление сбрасывается на $-\mathbf{g}$ (перезапуск) при обрезке PR+, при потере ортогональности градиентов и периодически каждые несколько итераций. Границы поиска параметров учитываются тем же штрафом; критерии останова — по норме градиента и по размеру шага.

6.8. Метод Брента — PRAXIS (Brent's method — PRAXIS), главных направленийбезградиентный

Объединяет метод поиска по направлениям (метод Пауэлла) с процедурой поиска вдоль главных осей (метод сопряжённых направлений с перестройкой базиса через SVD)^[8], обеспечивая быструю квадратичную сходимость. Мощный метод без производных, обычно наиболее эффективный на сложных и плохо обусловленных задачах.

Держит матрицу главных направлений и минимизирует вдоль них квадратичной аппроксимацией с дроблением интервала, делает квадратичный поиск по кривой второго порядка, а затем периодически перестраивает базис через SVD-разложение, выстраивая направления по собственным осям задачи; при плохой обусловленности включает случайные возмущения для выхода из «застревания».

Требует $n \ge 2$ (для одного параметра используется одномерный поиск Брента). Границы учитываются тем же штрафом; точность задаётся параметром $T_0$, длина шага подстраивается автоматически.

6.9. Метод QuickSearch — вращающиеся направлениябезградиентный

Безградиентный метод: вращающиеся направления с параболическим ускорением. Достаточно эффективный метод прямого поиска (вариант метода Розенброка) без производных. Держит набор из $n$ направлений и глобальный множитель шага.

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

Сходимость определяется по относительному размеру шага. Границы учитываются тем же штрафом.

6.10. Метод дифференциальной эволюции (jDE)глобальный

Алгоритм находит глобальный экстремум сложных математических функций, не требуя вычисления производных; метод популяционный (эволюционный)^[9]. Единственный глобальный метод в наборе данного проекта: остальные — локальные и сходятся в ближайший минимум. Целевая функция RSS может иметь несколько локальных минимумов (поверхность зависит не только от вида функции, но и от числовых значений исходных данных — коррелированные константы, недоопределённость); DE держит популяцию из NP векторов и умеет выходить из плохих бассейнов.

Механика (Storn — Price), каноническая схема DE/rand/1/bin — мутация:

$$ \mathbf{v} \;=\; \mathbf{x}_{r_1} + F\,(\mathbf{x}_{r_2} - \mathbf{x}_{r_3}) $$

(индексы $r_1$, $r_2$, $r_3$ различны и не совпадают с текущим), далее бинарный кроссовер с вероятностью CR и жадный отбор с элитизмом. Реализация усилена по трём осям:

Самоадаптация $F$ и CR (jDE, Brest 2006) — каждая особь несёт свои $F$, CR; они изредка перегенерируются и при удаче наследуются, снимая ручную настройку;
Мемётический полишинг — лучшая особь периодически доуточняется Гауссом — Ньютоном (RSS вблизи бассейна квадратична, GN полирует за единицы шагов), что резко ускоряет сходимость;
Финальное уточнение тем же Гауссом — Ньютоном перед расчётом доверительных интервалов — чтобы возвращаемая точка была настоящим минимумом (по ней строится гессиан).

Популяция засевается по всему боксу границ параметров; ГСЧ с фиксированным зерном (воспроизводимость).

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

6.11. Метод Ньютона второго порядка (демпфированный) (damped Newton)второго порядка

Использует не только градиент (первые производные), но и матрицу Гессе (вторые производные), что обеспечивает квадратичную сходимость — точность удваивается с каждой итерацией^[10]. Матрица Гессе описывает локальную кривизну многомерной функции и используется для точного определения характера критической точки как точки минимума.

Использовано демпфирование по Левенбергу и разложение Холецкого. Ньютоновский шаг:

$$ \delta \;=\; -\,\mathbf{H}^{-1}\mathbf{g}\,, \qquad \mathbf{g} = \nabla (RSS)\,, \quad \mathbf{H} = \nabla^{2} (RSS)\,, $$

где $\mathbf{H}$ — полный гессиан (в отличие от диагонального Левенберга — Марквардта и от $\mathbf{J}^{T}\mathbf{J}$ Гаусса — Ньютона), считается конечными разностями. У положительно определённого минимума это самый быстрый метод (квадратичная сходимость) — его можно применять напрямую на гладких задачах.

Вдали от минимума и в седловых точках (возможны при многоэкстремальности) $\mathbf{H}$ бывает знаконеопределена, и «чистый» Ньютон способен уйти в бесконечность. Страховка — демпфирование $(\mathbf{H} + \lambda \mathbf{I})$ с проверкой разложением Холецкого: успех означает положительную определённость (шаг — гарантированное направление спуска); при неудаче $\lambda$ повышается, пока матрица не станет положительно определённой (модифицированный Ньютон; $\lambda$ убывает у минимума $\rightarrow$ чистый Ньютон, растёт вдали $\rightarrow$ малый безопасный шаг, как у Марквардта). Вдоль направления — точный линейный поиск Брента.

Диагностика в точке решения. В найденной точке гессиан ещё раз проверяется Холецким: если он не положительно определён, выводится предупреждение (седло / овраг / неидентифицируемость — решение и доверительные интервалы ненадёжны).

Границы учитываются тем же штрафом; критерии останова — по норме градиента и по размеру шага.

Выбор метода и настройки

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

Флажок «Доуточнить методом Ньютона (второй порядок)». Завершает расчёт любым выбранным методом единообразно, что делает глубину минимума у разных методов сравнимой (иначе решения слегка различаются), и проверяет положительную определённость матрицы Гессе перед расчётом доверительных интервалов.

7.Доверительные интервалы

Решение задачи и доверительные интервалы параметров

Рис. 2. Решение той же задачи. Слева — набор параметров: отмечены E° и lg K(HAc), у каждого начальное значение и границы; справа — результат. Из заведомо сбитого старта (E° = 240 мВ, lg K = 4,0) Нелдер–Мид за 68 итераций пришёл к E° = 249,64 ± 2,6 мВ и lg K = 4,7653 ± 0,045 при истинных 250,0 и 4,76. Интервал здесь важнее самого числа: он говорит, насколько параметр вообще определён этими данными, а соседняя вкладка «Корреляции» — не определена ли вместо двух параметров лишь их комбинация.

Погрешности найденных параметров оцениваются по матрице Гессе целевой функции в точке минимума, вычисляемой численным дифференцированием (двухточечная схема по диагонали, четырёхточечная — для смешанных вторых производных).

Оценка дисперсии приближения и полуширины интервала:

$$ \sigma^{2} \;=\; \frac{F(\vec{\theta})_{\min}}{N - p}\,, \qquad \Delta\theta_i \;=\; t_{0{,}975;\,(N-p)} \sqrt{\sigma^{2}\,(\mathbf{H}^{-1})_{ii}}\,, $$

где $N$ — число точек, $p$ — число оптимизируемых параметров, $\mathbf{H}$ — гессиан, $t_{0{,}975}$ — квантиль распределения Стьюдента.

Вырожденность гессиана. При вырожденности гессиана интервал для соответствующего параметра не определяется. Это признак неидентифицируемости параметра по имеющимся данным — следует пересмотреть модель или набор уточняемых величин.

Корреляционная матрица — вкладка «Корреляции»

Доверительный интервал отвечает на вопрос, насколько плохо определён параметр сам по себе. Есть и второй вопрос, ответа на который по интервалам не видно вовсе: не определена ли вместо двух параметров лишь их комбинация. Оба интервала могут выглядеть приемлемо, а поверхность целевой функции быть длинным узким оврагом, вдоль которого решение гуляет свободно.

На этот вопрос отвечает корреляционная матрица, получаемая из того же обратного гессиана, что и интервалы:

$$ r_{ij} \;=\; \frac{(\mathbf{H}^{-1})_{ij}} {\sqrt{(\mathbf{H}^{-1})_{ii}\,(\mathbf{H}^{-1})_{jj}}} $$

Множитель $\sigma^2$ здесь сокращается, поэтому матрица достаётся даром — без единого дополнительного вычисления целевой функции. Она показана на вкладке «Корреляции» рядом с решением; ячейки с $|r| \ge 0.9$ подсвечены, с $|r| \ge 0.99$ — выделены красным. Прочерк означает, что величины не существует: дисперсия по этому параметру неположительна, то есть гессиан в найденной точке не положительно определён.

$\|r\|$	Что это значит	Что делать
< 0.9	параметры разделяются приемлемо	—
0.9 – 0.99	разделяются плохо: настоящий интервал каждого шире, чем следует из $\sigma$	расширить диапазон условий опыта; при публикации приводить корреляцию
≥ 0.99	практически неразделимы: данные определяют лишь комбинацию	закрепить один параметр независимыми данными либо не истолковывать значения по отдельности

В задачах этого комплекса высокая корреляция — норма, а не исключение. Типичные пары: $\lg K \leftrightarrow E^{\circ}$ (константа сдвигает расчётные активности почти равномерно по серии, а $E^{\circ}$ сдвигает всю кривую); параметры модели активности $\leftrightarrow$ обе эти величины. Именно поэтому переопределение $\mathring{a}$ и $b$ по отдельным ионам выключено по умолчанию, а перед его включением стоит посмотреть эту вкладку.

Матрица выгружается в Excel отдельным листом и копируется в буфер вместе с остальными таблицами окна.

Сравнение моделей активности по качеству подгонки

Вкладка «Сравнение моделей» подгоняет задачу заново каждой отмеченной моделью коэффициентов активности и раскладывает результаты рядом. Это второй, независимый от состава вопрос: какая модель лучше описывает эксперимент.

Сравнивать модели по сумме квадратов нельзя. Наборы параметров у моделей разные (см. 5.2 выше), а сумма квадратов монотонно падает с числом параметров — по ней «побеждала» бы всегда самая богатая модель. Поэтому мерой служит информационный критерий Акаике с поправкой на малую выборку, $\mathrm{AICc}$: меньше 2 единиц отставания — данные модели не различают, больше 10 — практически отвергают. Колонка RSS стоит рядом намеренно, и если порядок по ней расходится с порядком по AICc, программа говорит об этом отдельной строкой.

Подробный разбор — в справке «Расчёт равновесного состава», раздел 5.9; там же объяснено, как читать эту вкладку вместе со сравнением по составу.

8.Порядок и последовательность работы

Рис. 3. График невязок — главная проверка того, что модель верна по форме. Случайные знаки и отсутствие хода означают, что осталась только погрешность измерения; правильный ход невязок (дуга, ступень, рост к краю диапазона) указывает на пропущенную частицу, неверную модель коэффициентов активности или систематическую ошибку опыта, и никаким уточнением параметров не лечится. Ось X выбирается списком: номер раствора, любая аналитическая концентрация или ионная сила.

Ввод исходных данных

Шаг 1. Растворитель Для начала работы выбирается режим исследования «Потенциометрия», затем в меню «Файл» → «Новый» (Ctrl+N) для создания новой задачи, или «Открыть» (Ctrl+O), если требуется повторить расчёты с изменениями модели, заменить условия расчёта или на базе готовой модели создать новую. Это окно загружается во всех случаях.

Заполните поля «Название растворителя» (по умолчанию — вода), «Температура, К» (по умолчанию 298,15), «Диэлектрическая проницаемость» (по умолчанию для воды при этой температуре 78,303), «Плотность» (по умолчанию для воды при этой температуре 0,99707). По этим вводным автоматически рассчитываются параметры Дебая — Хюккеля для коэффициентов активности ионов в молярной шкале. При замене шкалы концентраций на моляльную в окне снизу данные будут пересчитаны с учётом плотности растворителя. Вводятся параметры уравнения Дебая — Хюккеля: параметр максимального сближения ионов $\mathring{a}$ и эмпирический параметр $b$.

При нажатии кнопки «Далее» — переход к продолжению ввода или смены информации. Кнопка «Назад» в первом окне не функционирует.

Шаг 2. Компоненты и матрица состава Укажите в окнах число компонентов и количество базисных частиц, после чего нажмите «Построить матрицу». Для ввода соответствующих названий (рекомендуется подробно) необходимо вызвать контекстное меню правой кнопкой мыши, щёлкнув по соответствующему полю заголовка. Далее ввести соответствующие элементы компонентной матрицы. Допускается ввод из буфера обмена.

Если требуется дополнить уже готовую матрицу строками или столбцами, необходимо изменить размерность и нажать «Построить матрицу» — при этом ранее введённая информация сохраняется, что удобно для формирования новых задач через готовые шаблоны.

Шаг 3. Стехиометрическая матрица Алгоритм ввода подписей и элементов матрицы аналогичный. Число строк, соответствующее числу частиц в растворе, вводится через окно ввода и пункт «Построить матрицу». Одновременно строится вектор зарядов частиц и вектор логарифмов констант равновесий образования соответствующих частиц.

Для частиц базиса элементы стехиометрической матрицы и вектора констант равновесий заполняются автоматически, так как равновесия базисных частиц формальны. Заполните все характеристики частиц и нажатием «Далее» перейдите к следующему окну.

Шаг 4. Концентрации компонентов Концентрации компонентов можно вводить вручную (для этого надо активировать «флажок» в поле окна) и через компонентную матрицу. Строки — растворы (по номерам), столбцы — компоненты (или частицы базиса, если включён прямой ввод). При большом числе растворов предусмотрена вертикальная прокрутка таблицы (колесо мыши / полоса справа).

Предусмотрена вставка данных в таблицу из буфера обмена (Ctrl+V): для этого необходимо выделить левую верхнюю ячейку и вставить таблицу целиком (при совпадении размерности) или начать с нужного столбца для частичной вставки.

На этом ввод общей информации для проведения вычислений закончен, и нажатием клавиши «Готово» производится переход к специфической для каждого раздела программного комплекса работе в зависимости от выбранного типа расчёта.

Предварительная оценка

При планировании эксперимента можно провести предварительный расчёт системы, перейдя в режим «Предварительная оценка». Это позволит, минуя процесс поиска решения оптимизацией параметров, рассмотреть численные характеристики: равновесные концентрации, модельные значения ЭДС гальванического элемента, значения коэффициентов активности и ионную силу растворов, а также сохранить рассчитанные результаты в файл в формате CSV.

Как и в режиме работы раздела «Равновесия», можно построить график зависимости равновесных концентраций (любых по выбору из всего числа частиц в системе), логарифма активности частиц базиса или коэффициентов активности в зависимости от выбранного пользователем аргумента: номер раствора, концентрация выбранного компонента, ионная сила.

Файл с введёнными исходными данными можно сохранить для повторного расчёта или в качестве шаблона в меню «Файл» (Ctrl+S), выбрав место его хранения и название.

Ячейка и ЭДС

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

Выберите в меню справа потенциалопределяющие ионы из общего списка частиц в растворе.

Концентрации компонентов, введённые или рассчитанные на предыдущем этапе, появятся в таблице автоматически. Введите экспериментальные значения ЭДС вручную или вставьте из буфера обмена. При отсутствии повторных измерений растворов веса автоматически приравнены 1; при вводе массивов повторных измерений будут рассчитаны значения весовой функции. Нажатие «Оптимизация» (F9) переводит в следующее окно для проведения процесса выбора и уточнения параметров равновесий.

Работа с окном оптимизации

В окне оптимизации в появившейся таблице, в которой представлены все возможные параметры, влияющие на функцию отклика, отметьте «галочкой» параметры, выбранные Вами для уточнения (в том числе нужные $\lg K$ — строки соответствуют частицам по стехиометрической матрице), укажите их начальные значения и, при необходимости, границы. Интервал возможного изменения параметров формируется автоматически как область в % от начального приближения, которое задаётся в окошке над таблицей. В случае необходимости изменения границ возможна ручная корректировка.

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

Нажмите «Оптимизировать». На правой вкладке отобразятся:

Вкладка	Содержание
Решение	найденные параметры $\vec{\theta}$ и доверительный интервал, представленный также нижней и верхней границей доверительной области
Точки	концентрации компонентов, экспериментальные значения $E_{\text{эксп}}$, веса $w_k$, рассчитанные значения ЭДС $E_{\text{расч}}$, невязки $\varepsilon_k$ по каждому раствору
Сходимость	график изменения $F(\vec{\theta})$ в зависимости от количества проведённых итераций
График невязок	информативно показывает качество минимизации и адекватность решения; аргумент зависимости выбирается по усмотрению в зависимости от цели анализа

Таблицы можно копировать в буфер обмена. Подробный анализ результатов полученного решения делается в меню правой вкладки.

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

Результаты можно экспортировать в Excel. Кнопка «Изменить частицы/равновесие» возвращает к стехиометрии для смены/добавления частиц.

9.Список литературы

[1]Nelder J. A., Mead R. A simplex method for function minimization. — Computer Journal, 1965, vol. 7, p. 308–313. Исходная работа о деформируемом многограннике: минимум ищется перестройкой симплекса из $n+1$ вершин — отражением, растяжением и сжатием, без единой производной. Метод §6.1 и умолчание приложения: он не требует ни гладкости, ни удачного начального приближения, а на овражной поверхности «$\lg K$ против $E^0$» это решает больше, чем скорость сходимости.
[2]Hooke R., Jeeves T. A. «Direct search» solution of numerical and statistical problems. — Journal of the ACM, 1961, 8 (2), p. 212–229. Работа, которая ввела само понятие прямого поиска: исследующий шаг по каждой координате плюс шаг по образцу вдоль накопленного смещения. Метод §6.2; собственный шаг по каждому параметру важен здесь потому, что $E^0 \approx 250$ мВ и $\lg K \approx 4$ — одного общего шага для них не существует.
[3]Levenberg K. A method for the solution of certain nonlinear problems in least squares. — Quart. Appl. Math., 1944, 2, p. 164–168. Демпфирование нормальных уравнений параметром $\lambda$: шаг непрерывно переходит от ньютоновского к наискорейшему спуску по мере ухудшения обусловленности. Метод §6.3 — тот самый приём, которым он и держится вдали от минимума. Там же по этому методу: Marquardt D. W. An algorithm for least squares estimation of nonlinear parameters. — SIAM J. Appl. Math., 1963, 11, p. 431–441 (правило изменения $\lambda$ по успеху шага — в этом виде метод и запрограммирован).
[4]Gratton S., Lawless A. S., Nichols N. K. Approximate Gauss-Newton methods for nonlinear least squares problems. — SIAM J. Optim., 2007, vol. 18, no. 1, p. 106–132. Разбор приближённого метода Гаусса — Ньютона и условий его сходимости. Метод §6.4 — единственный, которому нужен ЯВНЫЙ вектор невязок: из него строится якобиан, а $J^{\mathsf{T}}J$ передаёт связь параметров между собой, чего диагональное приближение §6.3 не делает.
[5]Nocedal J., Wright S. J. Numerical Optimization. — 2nd edition. — USA: Springer, 2006. — 664 p. Общее руководство по численной оптимизации; здесь — источник изложения BFGS (§6.5): обратный гессиан накапливается поправками ранга 2, а условие кривизны $y^{\mathsf{T}}s > 0$ решает, обновлять его или сбросить. Оттуда же взяты условия, которым отвечает линейный поиск вдоль направления. Там же по этому методу: Davidon W. C. Variable metric method for minimization. — SIAM Journal on Optimization, 1991, 1 (1), p. 1–17 (переиздание отчёта 1959 г. — работа, с которой началась переменная метрика и всё квазиньютоновское семейство); Avriel M. Nonlinear Programming: Analysis and Methods. — Dover Publishing, 2003.
[6]Powell M. J. D. An efficient method for finding the minimum of a function of several variables without calculating derivatives. — The Computer Journal, 1964, № 1 January (vol. 7), p. 155–162. Сопряжённые направления без производных: за внешнюю итерацию набор направлений обновляется составным шагом. Метод §6.6. Отсюда взята классическая форма проверки Флетчера — в приложении она восстановлена именно по первоисточнику, потому что без неё набор вырождается в координатные оси и поиск застревает в искривлённом овраге. Там же по этому методу: Powell M. J. D. Convergence properties of a class of minimization algorithms. — in: Nonlinear Programming 2, O. L. Mangasarian, R. R. Meyer, S. M. Robinson, eds. — Academic Press, 1977; Brent R. P. Section 7.3: Powell’s algorithm. — in: Algorithms for minimization without derivatives. — Englewood Cliffs, N. J.: Prentice-Hall, 1973. ISBN 0-486-41998-3 (линейный поиск Брента — общий для шести методов приложения).
[7]Polak E., Ribière G. Note sur la convergence de méthodes de directions conjuguées. — Revue Française d’Automatique, Informatique, Recherche Opérationnelle, 1969, 3 (1), p. 35–43. Формула $\beta$ Полака — Рибьера, сходящаяся заметно лучше классической Флетчера — Ривса. Метод §6.7; в приложении взят вариант PR+ с отсечением отрицательных $\beta$. Главное его достоинство здесь — память $O(n)$: обратный гессиан BFGS не хранится вовсе.
[8]Gegenfurtner K. R. PRAXIS: Brent’s algorithm for function minimization. — Behavior Research Methods, Instruments, & Computers, 1992, 24, p. 560–564. Описание PRAXIS в виде, пригодном для переноса в программу. Метод §6.8, самый тяжёлый в приложении: набор направлений периодически перестраивается через SVD и выстраивается по главным осям задачи — приём для плохо обусловленных случаев, где остальные безградиентные методы буксуют.
[9]Storn R., Price K. Differential Evolution — A Simple and Efficient Adaptive Scheme for Global Optimization over Continuous Spaces. — Technical Report TR-95-012, ICSI, March 1995. — 12 p. Дифференциальная эволюция — единственный ГЛОБАЛЬНЫЙ метод приложения (§6.10): популяция пробных точек, мутация разностью двух случайных особей. Берётся тогда, когда у задачи подозревают несколько минимумов; цена — полный расчёт равновесия на каждую пробу, то есть на один-два порядка дороже локальных методов. Там же по этому методу: Storn R., Price K. Differential Evolution — A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces. — Journal of Global Optimization, 1997, vol. 11, p. 341–359 (журнальный вариант); Price K., Storn R., Lampinen J. Differential Evolution: A Practical Approach to Global Optimization. — Springer, 2005. — 539 p. (самонастройка $F$ и $CR$, реализованная в приложении как jDE).
[10]Bonnans J. F., Gilbert J. C., Lemaréchal C., Sagastizábal C. A. Numerical optimization: Theoretical and practical aspects. — Universitext (Second revised ed. of translation of 1997 French ed.). — Berlin: Springer-Verlag, 2006. — 490 p. Источник изложения демпфированного метода Ньютона (§6.11): разложение Холецкого полного гессиана служит ОДНОВРЕМЕННО проверкой его положительной определённости, а неудача разложения поднимает демпфирование. Отсюда же — смысл диагностики: гессиан, не положительно определённый в найденной точке, означает, что доверительные интервалы недостоверны.