авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 10 |

«Ю.В.Холин КОЛИЧЕСТВЕННЫЙ ФИЗИКО- ХИМИЧЕСКИЙ АНАЛИЗ КОМПЛЕКСООБРАЗОВАНИЯ В РАСТВОРАХ И НА ПОВЕРХНОСТИ ХИМИЧЕСКИ МОДИФИЦИРОВАННЫХ КРЕМНЕЗЕМОВ: ...»

-- [ Страница 5 ] --

В некоторых программах использовали и другой метод нулевого порядка – метод деформируемого многогранника (симплекса) [250, 323, 324]. Симплексом называют правильный многогранник, образованный путем соединения отрезками прямых линий p + 1 точек (вершин многогранника), упорядоченных таким образом, чтобы для крите риальной функции выполнялись условия U1 U2... Up+1. При p = 2 симплекс – это треугольник, при p = 3 – пирамида. На каждой итерации текущий многогранник заме няется новым: точка xp+1 с наибольшим значением функции заменяется другой. Для p = 2 работу алгоритма легко представить наглядно (рис. 9.5).

Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности B J H G E A C Рис. 9.5. Работа симплекс-метода при p = 2.

Вначале случайным образом выбирают три стартовые точки – A, B и C, опреде ляющие исходный симплекс. В них оценивают критерий U. Пусть U(A) U(В) U(С).

Для создания нового симплекса выполняют следующие шаги:

1. Отбросить вершину A как имеющую наибольший отклик U(A).

2. Оценить U(E), где E – зеркальное отражение A:

3. Если U(E) U(C), оценить U(J), где J получают как еще один шаг вдоль вектора AE.

3.1. Если U(J) U(C), новый симплекс – BCJ. Перейти к шагу 6.

3.2. Если U(J) U(C), новый симплекс – BCE. Перейти к шагу 6.

4. Если U(B) U(E) U(C), новый симплекс – BCE. Перейти к шагу 6.

5. Если U(E) U(B), делается заключение, что многогранник слишком велик, и его надо сжать:

5.1. Если U(E) U(A), оценить U(G). Новый симплекс – BCG. Перейти к шагу 6.

5.2. Если U(A) U(E) U(B), оценить U(H). Новый симплекс – BCH. Перейти к шагу 6.

6. Начать построение нового симплекса.

В зависимости от характера поверхности минимизируемой функции по ходу ите раций происходит растяжение или сжатие симплекса. Итерации прекращают, когда вблизи минимума U многократное сжатие симплекса стягивает его в точку. Чтобы пре дотвратить получение не имеющих физического смысла отрицательных значений ис комых параметров, критериальной функции U, если один или несколько параметров становятся отрицательными, приписывают большое значение (скажем, 1050) [250].

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

Из методов нулевого порядка лучше всего зарекомендовал себя метод Хука Дживса [325]. Метод основан на следующей идее: в некоторой точке c пространства параметров находят, делая пробные шаги c±, направление движения, в котором ми нимизируемая функция уменьшается. Движение в этом направлении продолжают до той точки, в которой функция U или начинает расти, или перестает меняться. В этой точке направление движения выбирают заново. Метод обладает гарантированной схо димостью к локальному минимуму. Недостаток алгоритма Хука-Дживса, как и других методов нулевого порядка – чрезмерно низкая скорость сходимости.

Методы первого (градиентные) и высшего порядков используют информацию о величинах производных минимизируемой функции по искомым параметрам в некото рой точке c, чтобы выбрать направление движения из c в такую точку c+, что U(c+) U(c) (это не означает, конечно, что точка c+ обязательно лежит в на правлении градиента). Производные оценивают или как приращение конечных разно стей, или, что предпочтительнее, аналитически. На важную роль аналитического вы числения производных указывали, например, Т. Каден и А. Цубербюллер: «В програм мах не должно быть численного дифференцирования с его проблемами размера шага и потери значимых цифр. Преимущество аналитических производных – более высокая скорость и лучшая сходимость расчетов» [326].

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности Градиентные методы (антиградиента, скорейшего спуска, сопряженных градиен тов) [323, 327] основаны на том, что в произвольной точке c вектор градиента g функ ции U направлен в сторону возрастания U. В простейшем случае, минимизируя U, из меняют значения компонент вектора c в направлении антиградиента:

+ = c g ( c ), (9.76) где 0 1 – длина шага. Движение в пространстве параметров продолжают до тех пор, пока норма вектора + c не станет меньше заданной малой величины. Гради ентные методы могут сравнительно быстро прийти в окрестность минимума U. Однако если вблизи минимума функция имеет овражный характер, алгоритм, не останавлива ясь, блуждает вокруг этого минимума.

Наибольшее распространение в современных методах параметрической иденти фикации получили методы второго (квазивторого) порядка [323, 327, 328]. Важнейший из них – метод Ньютона. Его можно обосновать, рассмотрев квадратичную аппрокси мацию функции U(). Примем, что в окрестности некоторой точки c функция U не прерывно дифференцируема, и разложим ее в ряд Тейлора:

U ( ) = U ( c ) + g T ( c ) + ( c ) T H ( c ) +..., (9.77) где g – вектор градиента U, H – матрица Гессе производных U разме i j i ром pp, производные вычислены в точке c. Компоненты вектора g находят как N k U gi = = 2 wk k (9.78) i i k = или, в матричной форме, g = 2 JTW, (9.79) где W – диагональная матрица статистических весов wk размером NN, – вектор невя зок k, J – матрица Якоби производных Ak / i. Оставляя в разложении (9.77) три пер вых слагаемых, ищут вектор = - c, обращающий в минимум квадратичную форму Ф( ) = g T + ( ) T H ( ). (9.80) Поскольку = - с – стационарная точка функции Ф(), справедливо равенство H = -g = -2 JTW, (9.81) т.е.

p g ij j = gi, j = 1, 2,..., p. (9.82) i = Матрицу Гессе находят по формуле N wk ki kj + k i k j.

H ij = 2 (9.83) k = Выражение (9.81) эквивалентно формуле B = -JTW, (9.84) где B = JTWJ + F = (1/2) H, (9.85) F – матрица размером pp с элементами Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности N N 2 2 A wk k i j = wk k i j.

k k Fij = (9.86) k =1 k = В методе Ньютона для расчета вектора поправок решают систему нормальных уравнений = -B-1JTW. (9.87) Будь U() полиномом второй степени, движение от текущего значения с по вектору сразу, за один шаг привело бы в точку минимума. Для других функций введение попра вок + = c +, (9.88) где – длина шага, лишь уменьшает значение критерия U, но не приводит в ее мини мум. Полагая с := +, переходят к следующей итерации;

расчеты продолжают до схо димости.

Укажем некоторые важные характеристики метода Ньютона.

• Если матрица Гессе в точке минимума * положительно определена, а начальное приближение 0 лежит достаточно близко к *, то метод Ньютона сходится с квадра тичной скоростью1.

• Если квадратичная аппроксимация Ф() в окрестности точки c плохо воспроиз водит поведение функции U(), может наблюдаться расходимость итераций.

• Функция Ф() не имеет стационарных точек (в том числе конечной точки мини мума) в случае неположительно определенной матрицы Гессе H.

• При = 1 шаг вдоль Ньютоновского направления может привести к возрастанию U (хотя при этом достигается минимум аппроксимирующей функции Ф()!).

Алгоритмы, построенные на основе метода Ньютона, но компенсирующие его недостатки, называют модифицированными методами Ньютона. Из модифицирован ных методов Ньютона при параметрической идентификации моделей комплексообра зования чаще других использовали алгоритм Левенберга–Марквардта [323]. Если в ходе итераций в точке c матрица H стала неположительную определенной, ее моди фицируют:

~ H = H +I, (9.89) где – некоторое число, I – единичная матрица. Новую положительно определенную ~ матрицу H используют вместо Н в формуле (9.85). При алгоритм находит нап равление минимизации, параллельное вектору антиградиента -g.

В эпоху, когда ресурсы компьютеров были весьма ограниченными, метод Нью тона казался неэкономным, особенно при необходимости обрабатывать массивы дан ных, содержащие сотни или тысячи измерений. Основные затраты ресурсов ЭВМ свя заны с расчетом вторых производных, и метод Ньютона модифицировали таким обра зом, чтобы упростить (пусть и за счет точности) их оценивание. Основным квазинью тоновским алгоритмом является метод Гаусса-Ньютона [286, 323], внедренный впер вые в практику параметрической идентификации моделей комплексообразования Я. Ридбергом [329]. В методе Гаусса-Ньютона полагают, что в окрестности минимума невязки k и кривизна функции U() малы, и приравнивают нулю элементы матрицы F (см. уравнение (9.86)). Скорость сходимости метода близка к квадратичной, а затраты ресурсов ЭВМ значительно меньше, чем в методе Ньютона.

Еще одним квази-ньютоновским методом, нашедшем применение при расчете констант устойчивости, стал алгоритм Дэвидона–Флетчера–Пауэлла [250, 330]. В нем так же, как и в методе Гаусса-Ньютона, точную оценку матрицы Гессе заменяют при На каждой итерации число верных цифр в компонентах искомого вектора примерно удваивается (при ли нейной скорости сходимости каждая следующая правильная цифра требует в среднем трех итераций) [323].

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности ~ ближенной. Пусть на k-й итерации известна аппроксимация H k матрицы Гессе (исход ~ ное приближение H 0 обычно полагают равным единичной матрице). Переходя от k-й к (k+1)-й итерации, оценку матрицы Гессе уточняют:

~ ~ ~ H k +1 = H k + Bk, (9.90) где поправочная матрица 1 ~ ~ ~ H k ( k )( k ) T H k T T Bk = yk yk + ~ T ( k ) H k ( k ) yk ( k ), (9.91) T~ T + ( k ) H k ( k )µ k µ k yk = gk+1 - gk- приращение градиента на (k+1)-й итерации, 1 1 ~ µk = yk H k ( k ). (9.92) T~ T yk ( k ) ( k ) H k ( ) Метод обеспечивает сохранение положительной определенности матрицы Гессе при итерациях и не требует расчета вторых производных U. Впрочем, скорость i j сходимости метода к решению является линейной.

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

Для успешного расчета параметров равновесий необходимо использовать не один метод, а их комбинацию. Более того, предлагают минимизацию критериальной функ ции вести интерактивно, в диалоге «человек-компьютер» [331]. К сожалению, такой методологически правильной рекомендации весьма затруднительно следовать при мас совой обработке экспериментальных данных, в особенности если ее осуществляет спе циалист-химик, далекий от вычислительной проблематики.

9.3. Оценка адекватности модели Найдя минимум * критериального функционала U в пространстве параметров, исследователь переходит к оценке адекватности модели. Проверить адекватность мо дели – значит установить, описывает ли найденная модель экспериментальные данные с учетом погрешности их измерения.

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

9.3.1. Критерии 2 и F Это глобальные критерии. Если взвешенные невязки k = w1 / 2 k представляют k собой независимые случайные нормально распределенные величины с нулевым сред ним и единичной дисперсией, то такая случайная величина, как остаточная дисперсия N wk 2, f = N - z, s0 = (9.93) k f k = Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности распределена как /f с f степенями свободы и математическим ожиданием 1 [292, 332].

Если остаточная дисперсия чрезмерно велика, модель не адекватна эксперименту;

очень малая величина s 0 свидетельствует о переоценке уровня экспериментальных по грешностей. Модель считают адекватной, если выполняется двойное неравенство 2 (1 / 2) 2 = s0 f 2 ( / 2), (9.94) f f эксп где 2 (1 / 2), 2 ( / 2) – это 100(1-/2) и 100(/2)-процентные точки распределения f f 2 для f степеней свободы при заданном уровне значимости (обычно 0.10 или 0.05) [333]. Если приемлемы любые модели с s0 порядка единицы или меньше, проверяют справедливость неравенства 2 2 эксп = s0 f f ( ). (9.95) Соотношения (9.94), (9.95) точны только для моделей, линейных по параметрам.

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

Адекватность модели проверяют и с помощью F-критерия, сравнивая остаточную дисперсию s0 с независимой выборочной оценкой s2 дисперсии свойства А [266, 293].

Чтобы найти s2, необходимо выполнить повторные измерения в одной или нескольких точках плана эксперимента. Дисперсию s2 находят по результатам n параллельных из мерений свойства А как n (Ai A )2, w s2 = (9.96) n 1 i = где w – статистический вес A, рассчитанный так же, как и для других измеренных Ak.

Тогда случайная величина s F= 0 (9.97) s подчиняется F-распределению с числом степеней свободы ост = N - z и 0 = n - 1 [332, 333]. Если при заданном уровне значимости справедливо неравенство F Fтабл(;

ост;

0), (9.98) где Fтабл – 100-процентная точка F-распределения для числа степеней свободы ост и 0 [333], модель признают адекватной эксперименту. Если значение s2 принимают апри орно (при правильно назначенных статистических весах s2 = 1), то 0 =, а F-критерий совпадает с критерием 2.

9.3.2. Коэффициент асимметрии, эксцесс, средний остаток и среднее модулей остатков Зная величины взвешенных невязок k = w1 / 2 k, для распределения остатков k легко найти выборочные оценки коэффициента асимметрии ~~ ~ A = µ 3 (µ 2 ) 3 / 2, (9.99) ~ ~ где µ 2 и µ 3 – центральные моменты распределения остатков k (второй – дисперсия – и третий соответственно), эксцесса ~~ 2 = µ 4 (µ 2 ) 2 - 3, (9.100) ~ где µ 4 – четвертый центральный момент, среднего значения взвешенных остатков Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности N w1 / 2 k = (9.101) k N k = и среднего значения модулей взвешенных остатков N ~ w1 / 2 k.

= (9.102) k N k = Для расчета выборочных моментов используют формулы ~ µ = µ µ µ + 2(µ ) 3, (9.103) 3 3 21 ~ µ 4 = µ 4 µ µ 1 + 6 µ (µ 1 ) 2 3 (µ 1 ) 4, 3 2 (9.104) где N ( k ) j.

µ j = (9.105) N k = При нормальном распределении k с нулевым средним и единичной дисперсией мате матические ожидания указанных случайных величин составляют [333]:

~ ~ M ( A) = М(2) = М( ) = 0, M () = 2 / = 0.798. (9.106) ~ ~ Для проверки адекватности модели выборочные значения A, 2, и сравнивают с процентными точками распределения соответствующих статистик (табл. 9.4–9.6) [333].

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

~ Таблица 9.4. Процентные точки распределения статистики 100, % Число наблюдений 1 5 11 0.94 0.91 0. 16 0.92 0.89 0. 21 0.90 0.88 0. 26 0.89 0.87 0. 31 0.88 0.86 0. 36 0.88 0.86 0. 41 0.87 0.85 0. 46 0.87 0.85 0. 51 0.86 0.85 0. 71 0.85 0.84 0. 101 0.85 0.83 0. 201 0.83 0.82 0. 0.81 0.81 0. 9.3.3. Локальные критерии адекватности Замечено, что в ряде случаев модель, по глобальным критериям признанная адек ватной, в некоторых областях плана эксперимента неудовлетворительно описывает ре зультаты измерений. Более тонкую диагностику обеспечивают основанные на исследо вании взвешенных невязок k = w1 / 2 k локальные критерии. Проверяя адекватность с k ) ) их помощью, строят зависимость k от величин Ak. Примеры зависимостей k от Ak для адекватной и неадекватных моделей приведены на рис. 9.6–9.9.

Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности ~ Таблица 9.5. Процентные точки распределения статистики A 100, % Число наблюдений 1 25 1.06 0. 30 0.98 0. 35 0.92 0. 40 0.87 0. 50 0.79 0. 70 0.67 0. 100 0.57 0. 150 0.46 0. 200 0.40 0. 400 0.285 0. 1000 0.18 0. Таблица 9.6. Процентные точки распределения статистики 100, % Число наблюдений 1 50 1.92 1. 100 1.40 0. 200 0.98 0. 300 0.79 0. 500 0.60 0. 1000 0.41 0. k - Ak Рис. 9.6. Локальные критерии адекватности для адекватной модели.

out k - Ak Рис. 9.7. Визуализация резко выпадающего наблюдения.

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности k Ak Рис. 9.8. Систематический характер взвешенных невязок k для неполной модели.

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

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

Промахи стараются из обрабатываемых данных исключить. Однако даже при ис пользовании статистических критериев остается опасность неправомерно отбросить верное наблюдение или, наоборот, оставить в данных промах. Известно также [298], что далеко не всегда грубый промах проявляет себя в наибольшей невязке out. В связи с этим отбраковка выпадающих наблюдений на основе локальных критериев требует большой осторожности.

Более надежно обнаруживают локальные критерии неполноту модели (рис. 9.8).

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

Выявление таких областей – основа для пополнения модели новыми комплексом в ме тоде последовательной коррекции пробных моделей.

Визуализация остатков позволяет найти ошибки в назначении статистических ве сов [334]. Можно столкнуться со следующими ошибками: погрешность измерений за нижена (веса велики), завышена (веса малы) или неверно учтена зависимость погреш ностей от величин Ak. Рис. 9.9 иллюстрирует последний случай: веса wk были назна чены как 1 wk =, (9.107) 2 Ak s r где sr – относительная погрешность измерения Ak, одинаковая для всех Ak. На самом же деле при малых Ak погрешность s r s r, а при больших Ak s * s r.

* r k Ak - - Рис. 9.9. Локальные критерии адекватности при изменении sr с ростом Ak.

Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности 9.4. Проблемы единственности решений и избыточности моделей Пусть найдены значения параметров *, соответствующие минимуму минимизи руемого функционала, и построенная модель оказалась адекватной эксперименту. Те перь предстоит выяснить, не является ли модель избыточной, и не существуют ли дру гие адекватные модели.

В пространстве параметров в некоторой окрестности точки * существует замк нутая «область безразличия» * [335]. Все модели с параметрами, лежащими внутри этой области, также адекватны эксперименту. Построение области * при наличии развернутой информации о погрешностях измерений – рутинная операция прикладной математической статистики [266, 292]. Для поиска * используют ковариационную матрицу параметров D(*) (см. п. 9.5.1). Множественность решений, связанную с нали чием случайных погрешностей в экспериментальных данных, называют стохастиче ской [335].

Интереснее другой тип множественности решений – детерминированная неедин ственность [335], подразумевающая существование по крайней мере двух различных наборов параметров, одинаково хорошо (с учетом экспериментальных погрешностей) описывающих результаты измерений. Детерминированная неединственность возни кает, если минимизируемый критерий U() имеет в пространстве искомых параметров несколько изолированных минимумов (рис. 9.10). Обнаружить наличие нескольких ло кальных минимумов часто помогает повторение процесса минимизации, начиная с да леких друг от друга начальных оценок 0.

U * * 1 Рис. 9.10. Детерминированная неединственность модели при наличии нескольких локальных минимумов минимизируемого критерия.

При моделировании равновесий легко столкнуться с детерминированной неедин ственностью решений [273, 336–338]. Рассмотрим характерный пример из работ [273, 336], авторы которых по данным pH-метрического титрования построили две адекват ные модели комплексообразования в системе Cu2+-dl-лизин (HL). Каждая из моделей содержала четыре вида комплексов (табл. 9.7).

Предпочесть одну из моделей другой на основе данных рН-метрического экспе римента невозможно. Для дискриминации моделей необходимы дополнительные изме рения, связанные с переходом к другому диапазону общих концентраций компонентов [337] и/или измерению других свойств системы [336]. В рассматриваемом примере по могло привлечение нового экспериментального метода. Расчеты равновесного состава показывали, что в щелочной области для моделей 1 и 2 существенно различаются вели чины pCu (рис. 9.11) [273]. Следовательно, перспективным для дискриминации моде лей выглядит pCu-метрическое титрование. Его результаты позволили найти объеди ненную модель с пятью комплексами (табл. 9.7) [273].

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности Таблица 9.7. Модели комплексообразования Cu 2+ с dl-лизином Измеряемое Модель 1 Модель свойство lg 1 (Cu 2+ + L- = CuL +) lg 2 (Cu2+ + 2 L- = CuL2 ) pH lg 2 (Cu2+ + 2 L- = CuL2 ) lg 3 (Cu2+ + HL = CuHL2+) lg 3 (Cu2+ + HL = CuHL2+) lg 4 (Cu2+ +2 HL = Cu(HL)2 2+) lg 5 (Cu 2+ +L- + HL = CuL(HL)+) lg 4 (Cu2+ +2 HL = Cu(HL)2 2+) Umin = 24.8 Umin = 12. pCu объединенная модель:

lg 1 = 10.95 ;

lg 2 = 15.28 ;

lg 3 = 7.71 ;

lg 4 = 14.02;

lg 5 = 15. pCu pH Рис. 9.11. Кривая pCu-метрического титрования. Точки – эксперимент, 1 – расчет для модели 1, 2 – расчет для модели 2.

Расчеты равновесного состава для этой модели объяснили, почему данные pH метрии одинаково успешно описывают две модели: области максимального выхода комплексов CuL+ и CuLHL+, которыми отличаются модели, совпадают (рис. 9.12).

Кроме уточняющих экспериментов, возможна дискриминация адекватных моде лей при анализе соответствия каждой из них системе ранее полученных сведений о комплексообразовании [306, 339, 340].

Особые трудности в выявлении и преодолении детерминированной неединствен ности возникают в том случае, когда критериальная функция U в окрестности мини мума имеет овражный характер: существует по крайней мере одна линейная комбина ция параметров x = 11 + 2 2 +... + z z, (9.108) где i – числовые коэффициенты, для которой производная U / x = 0. (9.109) Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности [Li] pH Рис. 9.12. Зависимость равновесных концентраций комплексов ([Li ], ммоль/л) от рН для объединенной модели.

Как следствие, матрицы Якоби и Гессе являются матрицами неполного ранга, и не все искомые параметры поддаются определению. В терминах координационной хи мии это означает, что модель содержит избыточную химическую форму Lr (r – от re dundant, избыточный) с нулевым вкладом в измеряемое свойство A и материальный ба ланс [273, 284]. Стандартные средства линейной алгебры легко выявили бы и устрани ли избыточную форму Lr, будь ее вклад в измеряемое свойство и материальный баланс точно равен нулю. На практике же равенство (9.109) справедливо лишь приближенно, а соответствующие собственные числа матриц Якоби и Гессе приобретают малые, но не нулевые значения. В результате процесс минимизации теряет численную устойчивость.

Минимизация может завершиться и благополучно (итерации будут остановлены по ка кому-либо критерию, заложенному в алгоритм). Тогда реакции с участием избыточной формы Lr будет приписана некоторая константа равновесия Kr. Хотя фактически кри териальная функция U к значению Kr не чувствительна, потерявший бдительность ис следователь воспримет эту фиктивную константу как точно рассчитанную по результа там его измерений.

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

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

Начиная с работ Л.Г. Силлена, отбрасывали из модели те комплексы, константы устойчивости которых в процессе итераций становились меньше наперед заданной ма лой величины [261, 341, 342] или не превосходили своих средних квадратических от клонений [343]. При расходимости итерационного процесса Т. Каден и соавт. припи сывали одной из констант фиксированное значение, отказываясь от ее уточнения [344].

Л. Зекани и И. Надьпал отбрасывали из модели те комплексы, константы устойчивости которых на последней итерации менялись на одну логарифмическую единицу и более [262]. А.М. Евсеев и соавт. [324] и В.А. Бородин и соавт. [345] при обработке рН-мет рических данных устраняли из модели комплексы, равновесные концентрации которых во всех точках плана эксперимента не превосходили заданной пороговой величины.

Подобный подход при спектрофотометрическом исследовании равновесий использо вали авторы работ [250, 346], однако формы с малыми равновесными концентрациями Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности отбрасывали лишь в том случае, если рассчитанные молярные коэффициенты погло щения превосходили 105. Из методов борьбы с избыточностью моделей наиболее объ ективным можно считать метод, предложенный в работах [252, 316, 318]. Для выявле ния и устранения избыточных комплексов он использует свойства сингулярного раз ложения (SVD) [195, 203, 204] матрицы Якоби:

J = A / = U VT. (9.110) Алгоритм SVD заменяет рассчитанные поправки к логарифмам искомых констант lg j, j =1, 2, …, р, их линейными комбинациями p V ji lg j, xi = (9.111) j = где Vji – элементы ортогональной матрицы V, причем величины xi уже не коррелиро ванны, а их погрешности обратно пропорциональны сингулярным числам i – диаго нальным элементам матрицы. Принимая, что минимально возможная погрешность определения lg j составляет ~0.001, а незначимые lg j несут погрешность ~10-100, от ношение этих величин (~105) представляет собой предельную относительную погреш ность lg j [273]. Найдя отношения i / макс, приравнивают нулю те i, для которых это отношение не превосходит 10-5, и соответствующие поправки xi. В результате снижа ется размерность задачи, а неразрешимые (на данной итерации) линейные комбинации параметров не уточняются. Коэффициенты неразрешимых линейных комбинаций Vi, i = 1, 2, …, р, – это элементы столбцов матрицы V, номера которых соответствует номе рам малых сингулярных чисел [195]. Неразрешимая линейная комбинация логарифмов констант появляется тогда, когда в выражения для искомых констант входят концентра ции реагентов, не представленных в материальном балансе или измеряемом свойстве ни в одной экспериментальной точке. Для обнаружения избыточной химической формы достаточно просмотреть равновесные концентрации реагентов. В большинстве случаев и этого делать не приходится: уже по одним коэффициентам Vi удается выявить избы точную химическую форму. Недостаток метода состоит в том, что отбрасывание избы точного комплекса Lr происходит в процессе итераций, когда значения искомых пара метров еще далеки от точки минимума *. Возможно, что в окрестности минимума хи мическая форма Lr избыточной не является, и отбрасывание ее из модели в процессе расчетов было ошибкой.

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

В избыточности моделей проявляется математическая некорректность [79] задачи параметрической идентификации: «имеется много решений, совместимых с исходными экспериментальными данными, и ни одному из них нельзя отдать предпочтения: все они практически одинаково хорошо описывают исходные экспериментальные данные, если учитывать погрешность, с которой заданы последние» [347]. Принимая одно ре шение и отбрасывая остальные, мы руководствуемся не качеством описания экспери мента, а представлениями об искомых решениях. Рассмотренные выше алгоритм Ле венберга–Марквардта, SVD-разложение, а также гребневая регрессия [348] или регрес сия на ортогональных полиномах Чебышева [265] – это варианты Тихоновской -регу ляризации [79, 152, 212, 213].

Определенным подспорьем при выявлении детерминированной неединственно сти может служить разработанный в работах В.Г. Горского и Л.В. Быстрова [335] спо соб выявления нелинейных моделей неполного ранга. Он предполагает априорное, до выполнения эксперимента, исследование производных U /, сводящееся к выявле нию в матрице Якоби J линейно зависимых столбцов. При этом определяется ранг мо дели – максимальное число параметров, принципиально определимых в эксперименте рассматриваемого типа. Если исследуемая модель имеет неполный ранг, следующим шагом является нахождение таких нелинейных функций параметров i, которые и под Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности лежат оценке. Указанный способ анализа использует машинно-аналитические системы программирования (типа REDUCE), малая распространенность которых служила ос новным препятствием для внедрения этого подхода в практику моделирования. Крите риальные функции овражного типа, как указывают сами авторы [335], возникают чаще всего не как следствие неполноты ранга модели, а как результат неудачного планиро вания эксперимента. Тем не менее, априорный анализ ранга модели позволяет избе жать иллюзии, что выполнение повторных экспериментов, повышение их точности, пе реход к другим интервалам концентраций реагентов позволяет оценить принципиально неопределимые параметры. По этой причине нельзя не согласиться с рекомендацией «анализ ранга модели … проводить всегда, не считаясь со сложностью модели» [335].

9.5. Оценивание погрешностей параметров 9.5.1. Статистические оценки Точечные оценки параметров необходимо дополнить информацией об их по грешностях. Оценка погрешности параметров *, соответствующих минимуму крите рия U, основана на расчете ковариационной матрицы [286, 291, 292] [ ] D(*) = s 0 B ( * ), (9.112) где B – матрица, рассчитанная по формуле (9.85) при минимизации критерия U мето дом Ньютона, или ее приближенная оценка в квазиньютоновских методах. На диаго нали ковариационной матрицы содержатся оценки дисперсий компонентов вектора * – s 2 ( * ), а над и под диагональю – оценки ковариаций cov( *, * ) :

j i j s 2 (1 ) K cov(1, z ) cov(1, 2 ) s 2 ( 2 ) cov( 2, 1 ) L cov( 2, z ) D(*) =. (9.113) M M O M cov(, ) cov(, ) s 2 ( z ) L z1 z Матрица D(*) отличается от истинной ковариационной матрицы D(ист), т.к.

1) матрица B соответствует не истинным значениям искомых параметров ист, а их приближениям *;

2) при расчете ковариационной матрицы использована не истинная остаточная дисперсия s0, а ее выборочная оценка.

Если а) случайные величины s0 и * независимы (см. п. 9.2.2.1), б) случайная ве личина f s 0 (f = N - z – число степеней свободы) распределена как 2, в) распределение * ( ист ) i i * является многомерным нормальным, то случайная величина подчиня s ( i ) ется распределению Стьюдента [292], а интервал [ i t f, s ( i );

i + t f, s ( i )], (9.114) где tf, – множитель Стьюдента для f степеней свободы и доверительной вероятности, представляет собой индивидуальный 100%-ный интервал для (ист)i. В задачах нели нейного МНК свойства s 0 и *, использованные при построении интервала (9.114), имеют место лишь асимптотически при N. Чтобы сохранить для доверительных интервалов удобное выражение (9.114), вероятность считают заданной лишь прибли женно. Другой недостаток оценки (9.114) связан с тем, что она определяет доверитель ные интервалы каждого из параметров i по отдельности, не давая информации о веро ятности одновременного покрытия интервалами совокупности параметров. Вместо ин дивидуальных доверительных интервалов правильнее искать такую область R в про странстве параметров i, i = 1, 2, …, z, для которой вероятность Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности P(ист R) =. (9.115) Область R для приближенной доверительной вероятности 100 % представляют [266, 292] z-мерным эллипсоидом с центром в точке * и границами, задаваемыми ус ловием (* – )[D()]-1(* – )T = z F(1-, z, f), (9.116) где F(1-, z, f) – 100(1-)–процентная точка F-распределения для числа степеней свобо ды z и f. Доверительные эллипсоиды удобны при z = 2, но при z 2 теряют наглядность.

В этом случае предпочитают использовать совместные доверительные интервалы (ин тервалы Шеффе, Тьюки, Бонферони [292]). Удобен и прост метод Бонферони [267, 292], в котором индивидуальные 100-процентные доверительные интервалы параметров на ходят по формуле (9.114), принимая, что = 100 1. (9.117) z При высоких значениях z и интервалы Бонферони становятся чрезмерно боль шими. Например, при вычислении 95 %-ного индивидуального доверительного интер вала по формуле (9.114) при f = 10 множитель Стьюдента t = 2.23, а при построении интервалов Бонферони t = 3.17. Поэтому, оценивая интервалы Бонферони, рекоменду ется задавать совместную доверительную вероятность не выше 0.9.

При многомерном интервальном оценивании иногда пользуются привычной для химиков-аналитиков формулой t f, s ( i ) * ±, (9.118) i N взятой из приемов усреднения результатов параллельных измерений. Находимые по формуле (9.118) доверительные интервалы занижены (для логарифмов констант равно весия приводят [349] интервалы 0.001–0.005, тогда как, по оценке ИЮПАК, высокоточ ными признают значения lg, найденные с погрешностью 0.05-0.10 [283]). Выражение (9.118) является ошибочным, и применять его недопустимо.

Рассчитав ковариационные матрицы параметров, легко найти различные коэффи циенты корреляции [262]. Частные коэффициенты корреляции rij служат мерой взаим ного влияния параметров i и j при условии, что значения остальных параметров фик сированы:

Dij rij =. (9.119) ( Dii 1D 1 )1 / jj Общие коэффициенты корреляции sij измеряют степень коррелированности пара метров i и j, когда другие параметры рассматриваются как уточняемые:

Dij sij =. (9.120) ( Dii D jj ) Множественные коэффициенты корреляции Ri оценивают степень коррелирован ности параметра i со всеми остальными параметрами:

1/ Ri = 1. (9.121) Dii Dii Коэффициенты корреляции меняются в интервале [-1;

1]. При нулевых значениях параметры независимы, при значениях ±1 модель переопределена (существует линей ная зависимость между параметрами). Считают, что отбрасывание из модели одного из комплексов необходимо, если общий коэффициент корреляции превышает пороговое значение 0.95-0.98.

Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности 9.5.2. Нестатистический подход Неудовлетворенность тем, что статистические методы расчета параметров и оце нивания их доверительных интервалов целиком базируются на априорных гипотезах о статистических свойствах моделируемых объектов, побудила искать альтернативу.

Г.А. Коковин, В.А. Титов и соавт. обосновали [350, 351] достоинства предложенного Л.В. Канторовичем [352] подхода для математических задач химической термодина мики. Этот подход свободен от гипотез о распределении погрешностей эксперимен тальных данных, что позволило назвать его нестатистическим.

Пусть экспериментально изучена зависимость некоторого свойства A от контро лируемых переменных xj;

модель системы ~ A = f (xj;

) (9.122) задана с точностью до параметров ;

каждая экспериментальная величина измерена с погрешностью, максимальная величина которой известна:

э ист э Ak Ak ( Ak ) Ak + Ak, k = 1, 2,..., N, (9.123) ( x э ) k (x j ) k ( x ист ) k ( x э ) k + (x j ) k, k = 1, 2,..., N, (9.124) j j j где индексом «э» отмечены экспериментальные значения, а интервалы Ak, (xj)k – мак симальные погрешности Ak и xj в k-ой экспериментальной точке. С учетом функцио нальной зависимости (9.122) решают относительно искомых параметров систему нера венств (9.123), (9.124). Решение этой системы представляет собой область допустимых значений (ОДЗ) параметров. Систему неравенств можно дополнить условиями, вытека ющими из физического смысла задачи (например, i 0) и ограничивающими ОДЗ. К достоинствам метода относят следующие [350]: результаты надежны с метрологичес кой точки зрения, поскольку Ak и (xj)k принимают во внимание как случайные, так и систематические погрешности;

максимально возможным образом используется экспе риментальная информация об объекте исследования;

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

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

Пересмотр гипотезы о погрешностях в сторону их уменьшения снижает размер ОДЗ.

Зачастую точечные оценки, соответствующие минимуму критериальной функции U, лежат вне ОДЗ (рис. 9.13).

2 * Рис. 9.13. ОДЗ для двух варьируемых параметров. 1 – высокая оценка уровня экспериментальных погрешностей;

2 – низкая оценка;

3 – точечная оценка, найденная как решение задачи МНК.

К сожалению, до сих пор не решены проблемы представления ОДЗ и хранения информации для задач большой размерности;

не исключено сильное влияние на ОДЗ даже единичных грубых промахов;

серьезные вычислительные трудности возникают для многомерных нелинейных по параметрам задач [350]. Эти обстоятельства не по Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности зволяют рассматривать нестатистический подход как полноценную замену методу мак симума правдоподобия. Он полезен в той мере, в которой помогает более критично оценивать результаты статистического подхода.

9.6. Алгоритм CLINP 1.0 обработки данных многоволновой спектрофотометрии 9.6.1. Основы расчета параметров и их ковариационной матрицы Наш алгоритм CLINP 1.0 [252, 253] можно считать характерным представителем семейства алгоритмов нелинейного МНК, прошедших многолетнюю апробацию при моделировании комплексообразования в сложных системах. Опишем кратко этот алго ритм, обращая особое внимание на те его особенности, которые обеспечивают числен ную устойчивость и скорость расчетов. Алгоритм предложен для обработки данных многоволновой спектрофотометрии, однако естественным образом распространяется и на методы потенциометрии, растворимости, распределения, ядерной магнитной релак сации и др.

Примем, что • известны число химических форм в растворе S и их стехиометрический состав;

• реакции в системе записаны в канонической форме и по условиям смешивания реагентов точно известны общие концентрации всех Y компонентов в N растворах;

• условия равновесия выражены уравнениями ЗДМ (1.5) с использованием концен трационных констант;

• измерено светопоглощение N растворов в условиях равновесия при длинах волн;

• для всех растворов и длин волн справедливы уравнения Фирордта [243] S Alk = li [Li ]k + lk, (9.125) i = где l – номер аналитической позиции (длины волны);

k – номер раствора;

Alk – светопоглощение k-го раствора при длине волны l, деленное на длину поглощающего слоя;

li – молярный коэффициент поглощения реагента Li при длине волны l;

[Li]k – равновесная концентрация Li в k-м растворе, lk – случайная погрешность измерения Alk. Величины Alk образуют матрицу A размером N, li – E размером S, [Li]n – C размером SN. Предположим также, что погрешности li – независимые случайные нор мально распределенные величины с нулевым средним и дисперсией 2. В матричной записи уравнения (9.125) имеют вид A = EC +. (9.126) Расчету подлежат: неизвестные константы устойчивости i, i = 1, 2, …, p;

p S, и факторы интенсивности реагентов, записанных в стехиометрической матрице пер выми. Общее число искомых параметров z = p + N. (9.127) Разобьем матрицы Е и С на блоки, соответствующие реагентам с известными и неизвестными факторам интенсивности:

Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности S E= E1 E N C C= p C S где E1 – блок размерностью с неизвестными факторами интенсивности li, Е2 – блок с известными l. Матрицу A представим в виде A = E1C1 + E2C2 +, (9.128) + A E1 E2 C = x C Параметры ln i, li определяем методом наименьших квадратов, минимизируя критериальную функцию ( )( ) N U = lk = tr T = tr T, (9.129) =1 k = где невязки lk = Alk Alk, – матрица невязок размером N, tr – оператор следа.

Светопоглощения Alk линейно зависят от li и нелинейно – через [Li] – от ln i.

Используя подход Л.Г. Силлена (см. п. 9.2.2.4) [250, 353], параметры li оцениваем ли нейным МНК, не требующим начальных приближений и сходящимся за один шаг. За дав некоторые оценки c, можно вычислить матрицы C и E1:

)( )1.

( E1 = AT E 2 C C1 C1C T T (9.130) Оценки E1 необходимы при итеративном уточнении ln c. Критерий U миними зируем методом Гаусса-Ньютона, в котором для псевдообращения матриц использо вано сингулярное разложение. Приближенные оценки (ln i)c, i =1, 2, …, p, итеративно уточняем, вводя к ним поправки i:

ln i + = (ln i ) c i, i = 1, 2,..., p, (9.131) где – длина шага. Поправки i находим, решая линейным МНК систему уравнений p ( lk ln i ) i = lk, l = 1, 2,..., ;

k = 1, 2,..., N. (9.132) i = Производные lk ln i составляют трехмерную матрицу D, D – оператор дифференцирования / ln. Матрицу D подвергаем SVD-разложению. Если при этом находим такие сингулярные числа i, для которых i / макс 10-5, то выявляем и отбрасываем неразрешимые на данной итерации линейные комбинации параметров Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности ln i, к варьированию которых критериальная функция U не чувствительна (см. п. 9.4).

Тем самым обеспечивается численная устойчивость расчетов и устранение избыточных параметров.

Длину шага в (9.131) ограничиваем:

p = min1;

15 i.

(9.133) i = Если при новых оценках ln i, полученных по формуле (9.131), критерий U убы вает, принимаем их. В противном случае дробим шаг пополам и снова применяем фор мулу (9.131). Дробление шага повторяем до тех пор, пока не добьемся убывания U.

Уточнив ln i, заново рассчитываем матрицы C, E1, и поправки i. Итерации p i2 превышает 10-5.

продолжаем, пока сумма квадратов поправок i = Поскольку матрица C в процессе минимизации вычисляется многократно, для расчета равновесного состава требуются очень быстрые алгоритмы. Из множества ис пытанных наиболее быстрым и надежным мы сочли проекционный метод А.А. Бугаевского [354].

Производные lk ln i вычисляем аналитически:

D = D (E1C1 + E 2 C 2 A) = (DE1 )C1 + E1 (DC1 ) + E1 (DC 2 ). (9.134) Производные DE1 получаем при дифференцировании выражения (9.130):

)( ) ) ) ( ( ( T T 1 T T DE1 = ( A E 2 C 2 ) DC1 C1C1 E 2 (DC 2 )C1 C1C T T + C1 D C1C1 (9.135) где )[ )](C1C1T )1.

) ( ( ( T 1 T (DC1 )C1 + C1 DC T T = C1C1 (9.136) D C1C Окончательно выражение для производных принимает вид:

D = DA = ( DC1T )µ T + E ( DC ) (1 µ C1 ), (9.137) где µ = C1T(C1C1T)-1. (9.138) Необходимые производные DC1 и DC2 получены впервые А.А. Бугаевским и со авт. [312, 313]:

Y S Y [ Li ] / ln h = [ Li ] ih - ij dh du Puj [ Ld ], (9.139) j =1 d =1 u = где – символ Кронекера, Puj – элементы матрицы буферности [275], || Puj || = || [Bu] / tj || = || tj / [Bu] ||-1. (9.140) По окончании уточнения ln i рассчитываем, следуя [317, 355], ковариационную матрицу параметров ln i:

( ) [( ) ] cov ln = s0 D T lki * (D )lki, l = 1, 2,..., ;

k = 1, 2,..., N ;

i = 1, 2,..., p, (9.141) s0 = U / f, f=N-z (9.142) В формуле (9.141) суммирование выполняется по повторяющимся индексам l и k.

Ковариационная матрица совокупности параметров ln i и li имеет блочную форму Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности cov (ln ) cov cov ( ) =, (9.143) cov T cov 3 где матрицы cov2, cov3 построены из блоков, соответствующих длинам волн l: cov, cov 3, l = 1, 2, …,. Каждый из блоков cov 2 представляет собой матрицу размером, содержащую дисперсии–ковариации li для длины волны l:

)1 + (cov 3 )T (D ).

( cov = s0 C1C 2 T (9.144) Блок cov 3 размером p содержит ковариации ln i и li для длины волны l:

(cov 3 )T = (D ) cov(ln ), (9.145) 9.6.2. Снижение размерности задачи по числу аналитических позиций Для многоволновой спектрофотометрии характерны значения от 5–10 до 102– 10. При большом числе длин волн расчетные алгоритмы должны содержать средства, снижающие размерность задачи без потери существенной информации. Одно из таких средств – факторный анализ [155, 258], основанный на сингулярном разложении мат рицы A:

A = U V T. (9.146) При высоких и ранг матрицы А (r) не превышает числа светопоглощающих форм. Ранг A равен числу сингулярных чисел (jA), значимо отличающихся от 0. Такими признавали те числа (jA), для которых (jA) / ( A) 0.01.

макс Почти вся информация, содержащаяся в матрице A, теперь сохраняется в заштри хованных частях матриц U и V:

r r N r rT V N U Это позволяет заменить A меньшей по размеру матрицей – произведением за штрихованных частей VT. К ней переходим, умножая соответствующие матрицы на заштрихованную часть UT. Аналогичному преобразованию подвергаем E, E1, E2, и D.

9.6.3. Схема алгоритма CLINP 1. 1. Выполнить сингулярное разложение матрицы A, найти ее ранг r и преобразо ванные матрицы VT = UTA и UTE2.

2. Задать начальные приближения (ln i)0, i =1, 2, …, p.

3. Для текущих значений (ln i)c рассчитать матрицы C и UTE1.

4. Линейным методом наименьших квадратов решить систему уравнений (9.132) и рассчитать поправки i.

5. Ввести поправки i к текущим оценкам (ln i)c.

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности p i2 10 5, перейти к шагу 3, в противном случае – к шагу 7.

6. Если i = 7. По обратному U-преобразованию найти матрицы E1, и D. Вычислить ко вариационную матрицу cov().

9.6.4. Пример использования алгоритма CLINP 1.0. Равновесия в системе Cu(II) – этилендиамин – оксалат-ион На основе описанного алгоритма разработана программа CLINP 1.0 [253]. Приве дем пример ее использования при моделировании равновесий комплексообразования в растворах. В монографии [244] подробно исследованы равновесия в растворе, содержа щем Cu(II), этилендиамин (En) и щавелевую кислоту (Н2Ох). При шести длинах волн и 21 значении рН измерены светопоглощения растворов с общими концентрациями t(Cu2+) = 0.01 моль/л, t(En) = 0.1 моль/л, t(H2Ox) = 0.01 моль/л. Авторы [244] анализиро вали равновесия в этой системе как с применением вспомогательных функций и гра фических методов, так и с помощью компьютерного моделирования по программе DALSFEC [244, 248]. Оказалось, что поглощают свет три частицы: CuEn22+, CuEnOx и CuOx22-, а концентрация ионов Сu2+ пренебрежимо мала. Рассчитаны K1 и K2 – кон станты равновесия реакций K CuEn 2+ + Ox 2 = CuEnOx + En, (9.147) K CuEnOx + Ox 2 = CuOx 2 + En (9.148) и молярные коэффициенты поглощения комплексов (li) [244].

По программе CLINP 1.0 рассчитывали lg K1, lg K2, lg 2( Cu 2 + + 2 En = CuEn 2 + ) и li. Результаты расчетов по программам CLINP 1.0 и DALSFEC приведены в табл. 9.8.

Таблица 9.8. Результаты моделирования равновесий в системе Cu 2+ – En – Ox2-.

Верхняя строка – расчет по программе DALSFEC, нижняя – по программе CLINP 1. -lg K1 -lg K U s2 * * 4.79(0.02) 5.96(0.02) 1.2110 -3 - 4.78(0.02) 5.95(0.02) 1. Коэффициенты молярного поглощения комплексов Длина волны, нм CuEn 2+ CuOx CuEnOx 2 480 26.9 (0.2) 2.4 (0.3) 0.7 (0.2) 26.9 (0.1) 2.4 (0.2) 0.7 (0.2) 600 46.0 (0.2) 47.1 (0.35) 11.5 (0.2) 46.0 (0.1) 47.2 (0.2) 11.5 (0.1) 640 28.5 (0.2) 46.6 (0.35) 23.2 (0.2) 28.5 (0.1) 47.7 (0.2) 23.2 (0.1) 646 26.2 (0.2) 46.6 (0.35) 24.8 (0.2) 26.2 (0.1) 46.6 (0.2) 24.8 (0.1) 700 11.5 (0.2) 33.1 (0.3) 34.0 (0.2) 11.5(0.1) 33.1 (0.2) 34.1 (0.1) 800 2.5 (0.2) 11.8 (0.3) 24.85 (0.2) 2.5 (0.1) 11.8 (0.2) 24.9 (0.1) Глава 9. Модель химических реакций: опыт расчета параметров и проверки адекватности [Li] 10 1 5.0 5.5 6.0 6.5 7. pH Рис. 9.14. Рассчитанные равновесные концентрации (ммоль/л) комплексов CuEn 2 + (1), CuEnOx (2) и CuOx 2 (3) в исследуемых растворах.

2 Рассчитанная по программе CLINP 1.0 остаточная сумма квадратов U = 1.2110- при числе степеней свободы f = 106. Принимая оценку среднего квадратического отклонения светопоглощений A = 0.003 и уровень значимости = 0.025, приходим к выводу об адекватности модели по критерию 2:

2 = 1.21 10 3 /(3 10 3 ) 2 = 134.4 106 (0.025) = 136.4.

эксп (9.149) Рассчитанные равновесные концентрации комплексов (рис. 9.14.) объясняют, почему удалось надежно рассчитать все li: каждый из комплексов имеет широкую область до минирования. Параметр 2 не поддается оценке по имеющимся данным, т.к. химиче ская форма Cu2+ является избыточной (концентрация [Cu2+] не превышает 0.1% t(Cu2+) ни в одной экспериментальной точке).

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


Проанализировав результаты, достигнутые в практике параметриче ской идентификации моделей, можно заключить:

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

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

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности 3. К настоящему времени основными методами минимизации критериаль ного функционала U являются модификации метода Ньютона и метод Гаусса-Ньютона. Алгоритмы, ставшие основой распространенных про грамм, лишь в незначительной степени учитывают результаты, получен ные специалистами по нелинейной оптимизации в последние десятиле тия. Как результат, скорость и сам факт сходимости алгоритмов к реше нию в значительной степени зависят от выбора начальных приближений параметров и эвристических соображений;

не решена до конца проблема выявления и устранения детерминированной неединственности решений.

4. Оценка адекватности моделей по глобальным критериям 2, F и т.п. не достаточна. Необходимо анализировать локальные критерии адекватно сти и внедрять новые методы диагностики моделей.

5. Все выводы о доверительных интервалах и коррелированности парамет ров, основанные на расчете ковариационных матриц, являются прибли женными.

ГЛАВА 10. МОДЕЛЬ ХИМИЧЕСКИХ РЕАКЦИЙ:

ОЦЕНИВАНИЕ ПАРАМЕТРОВ ПРИ НАРУШЕНИИ СТАНДАРТНЫХ ПРЕДПОСЫЛОК МНК. М-ОЦЕНКИ ХЬЮБЕРА 10.1. Нарушения предпосылок МНК и меры по их компенсации В методе максимума правдоподобия выбор критериальной функции в виде N U = 2, (10.1) k k = где взвешенные невязки k = w1 / 2 k, k = Ak - Ak, основан на гипотезах о точном зада k нии всех переменных, контролирующих условия эксперимента, и о нормальном рас пределении (погрешностей измерения свойства A).

Первая гипотеза справедлива лишь приближенно, поскольку все эксперимен тальные величины, включая и те, которые задают условия эксперимента, подвержены случайным погрешностям. При нормальном распределении всех погрешностей оценки параметров с оптимальными свойствами (см. п. 9.2.2.1) следует находить с помощью конфлюэнтного анализа (регрессионного анализа с учетом ошибок в контролируемых переменных) [293]. В конфлюэнтном анализе вместо (9.47) минимизируют функционал N ~ ( U = wk k, (10.2) k = ( ) где отклонение k учитывает как невязку k = Ak Ak, так и ошибки в векторе конт ролируемых переменных Х [294]:

1 ( k = k + 2 d, (10.3) 2, k Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера где индексы и пробегают все значения, соответствующие контролируемым пере менным;

d – элементы матрицы погрешностей контролируемых переменных 2 A 2 =. (10.4) X X Приходится корректировать и значения статистических весов:

~ = s 2 + 2 d k 1 wk, (10.5), k где s k – оценка дисперсии Ak, A A 1 =. (10.6) X X Основная сложность в применении конфлюэнтного анализа – отсутствие досто верной информации о погрешностях контролируемых переменных. Немногочисленные примеры использования конфлюэнтного анализа [294–296] показывают, что конфлю энтные поправки незначительно (в пределах погрешности определения) влияют на оценки i. Это дает основание не пользоваться конфлюэнтным анализом, если случай ные погрешности контролируемых переменных невелики.

Гораздо сложнее проблемы, возникающие при нарушении гипотезы о нормаль ном распределении погрешностей. В этом случае МНК может дать абсолютно невер ные результаты. Рассмотрим простой модельный пример, подобранный П. Хьюбером [298]. В шесть откликов, соответствующих линейной модели y = 2 x, (10.7) Хьюбер внес погрешности: в y1 – y5 – нормально распределенные с нулевым средним и стандартным отклонением = 0.6, а в y6 – большую ошибку 12 (табл. 10.1).

Таблица 10.1. Результаты регрессионного анализа при распределении погрешностей, отличном от нормального (пример П. Хьюбера) No Модель (10.8) Модель (10.10) xi yi ) ) i i yi yi 1 -4 2.48 0.39 -2.09 2.23 -0. 2 -3 0.73 0.31 -0.42 0.99 0. 3 -2 -0.04 0.23 0.27 -0.09 0. 4 -1 -1.44 0.15 1.59 -1.00 0. 5 0 -1.32 0.07 1.39 -1.74 -0. 6 10 0 -0.75 -0.75 0.01 0. Расчет линейным МНК коэффициентов уравнения y = a0 + a1x (10.8) дал оценки a0 = 0.07 и a1 = -0.08, весьма далекие от верных значений. При этом оста точная дисперсия 16 i = 2.41, s0 = (10.9) 4 i = ) где i = y i y i, хотя и существенно превышает 2 = 0.36, пугающе большой не выгля дит. Исследование локальных критериев адекватности i указывает на первую точку как на возможный грубый промах (невязка |1| = 2.09 наибольшая), т.е. грубый промах не проявил себя в наибольшей невязке.

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности y 4 - - - - - - -4 -2 0 2 4 6 8 x Рис. 10.1. Результаты регрессионного анализа для примера П. Хьюбера.

1 – модельная функция, точки – данные с внесенными погрешностями, 2 – уравнение линейной регрессии (10.8), 3 – полином второй степени (10.10).

При внимательном взгляде на данные (рис. 10.1) возникают обоснованные со мнения в пригодности линейной модели (10.8). Полиномиальная модель y = a 0 + a1 x + a 2 x 2 (10.10) 2 приводит к МНК-оценкам a0 = -1.74;

a1 = -0.66;

a2 = 0.08 и s 0 = 0.17. Значение s 0 пре красно согласуется с дисперсией 2, и полиномиальная модель (10.10) будет принята.

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

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

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

Метод максимума правдоподобия приводит к различным критериальным функ циям в зависимости от распределения экспериментальных погрешностей. Так, при рав номерном распределении оценки максимума правдоподобия соответствуют мини муму максимальной невязки [356, 357] S = max k. (10.14) 1 k N Если распределение подчиняется закону Стьюдента, следует минимизировать крите рий [358] ( ) N Z = m k ln 1 + m k 2 v 1 s k 2, (10.13) k k k = где mk – число повторных измерений Ak, vk = mk - 1, s k – выборочная оценка дисперсии Ak. Наконец, при распределении, соответствующем закону Лапласа с плотностью ве роятности exp( ) p ( ) = (10.11) где – параметр распределения, критериальная функция имеет вид [299, 300, 359] Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера N | k |.

L1 = (10.12) k = Критерий S не годится, т.к. он еще более чувствителен к засорению данных про махами, нежели U. Оценки, соответствующие минимуму критериев L1 и Z, напротив, значительно более устойчивы к промахам, чем МНК-оценки. Вместе с тем, критерий Z не хотелось бы использовать по той причине, что закон Стьюдента не является пре дельным, и трудно представить экспериментальные условия, формирующие Стьюден товское распределение погрешностей. Напротив, как показали Мудров и Кушко [299] (см. п. 9.2.2.1), критерий L1 статистически обоснован так же хорошо, как и критерий U, т.к. закон распределения Лапласа, как и закон Гаусса, является предельным. Обосно ванно предпочесть один критерий другому в отсутствие информации о распределении экспериментальных погрешностей невозможно. Это тем более печально, что миними зация различных критериев может приводить к значительно отличающимся оценкам параметров [331, 357].

Решение проблемы следует искать, задавая минимизируемый критерий так, чтобы он: а) давал оценки, устойчивые к отклонению распределения погрешностей от нормального, и б) объединял два статистически обоснованных критерия – U и L1.

Статистические методы, устойчивые к нарушению стандартных предпосылок МНК, называют робастными [291, 298, 360, 361]. Робастные методы оценивания не только повышают устойчивость оценок к засорению данных промахами, но, как мы в этом убедимся, смягчают (по сравнению с МНК) требования к назначению статистиче ских весов и выбору той или иной характеристики равновесной системы в качестве ве личины, аппроксимируемой моделью. Поскольку при робастной параметрической идентификации приходится варьировать гипотезы о статистических свойствах экспери ментальных данных, есть все основания рассматривать робастное оценивание как раз дел теории анализа данных.

10.2. М-оценивание как робастный метод параметрической идентификации Показано, что оценки МНК теряют свои оптимальные свойства при распределе нии погрешностей с хвостами более длинными, чем у нормального распределения;

при более коротких хвостах МНК-оценки удовлетворительны [286]. Следовательно, конструируемый робастный критерий должен обладать свойствами взвешенной суммы квадратов при распределениях погрешностей с короткими хвостами и быть устойчи вым к распределениям с длинными хвостами. Количественно длину хвостов можно ха рактеризовать эксцессом 2 (9.100). Для нормального распределения 2 = 0, для распре делений с длинными хвостами 2 0 (например, для распределения Лапласа 2 = 3, для распределения Стьюдента 2 =, где f – число степеней свободы, для распределе f ния 2 эксцесс 2 = [332, 362]), для распределений с короткими хвостами эксцесс f отрицателен (так, для равномерного распределения 2 = -1.2 [332, 362]).


Строя робастный критерий методом максимума правдоподобия, применяют «мо дель грубых промахов» [293]. Пусть, свойство А измеряется со случайной погрешно стью, плотность распределения которой р() = [(100 - )(0, Гаусс) + h()] / 100, (10.18) где (0, Гаусс) – плотность стандартного нормального распределения (среднее значение 0, стандартное отклонение Гаусс = 1), h() – плотность вероятности грубых промахов (плотность распределения с длинными хвостами), – интенсивность промахов, %.

Примером модели грубых промахов служит известное распределение Тьюки [362] Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности р() = [(100 - )(0;

Гаусс) + (0;

bГаусс)] / 100, b 1. (10.19) В уравнении (10.19) h() = (0;

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

Гаусс).

П. Хьюбером показано [298, 360], что для модели грубых промахов оценка мак симального правдоподобия соответствует минимуму критериального функционала N M () = ( k ;

) = min (10.20) k = или решению уравнения N ( k ;

) = 0, (10.21) k = где (k;

) – неотрицательная, симметричная и дважды дифференцируемая функция потерь, для которой (0;

) = 0, ( k, ) = ( k )( k, ), k – взвешенные невязки. За дачи (10.20) и (10.21) эквивалентны, если функция выпуклая. Оценки, соответствую щие минимуму критерия (10.20), называют М-оценками. М-оценки несмещенные, сос тоятельные и эффективные, если функция (;

) растет медленнее, чем 2. М-оценки устойчивее МНК-оценок к нарушениям гипотезы о нормальном распределении экспе риментальных погрешностей [293, 363].

Предполагая, что функция плотности вероятности грубых промахов h() симмет рична относительно нуля, Хьюбер потребовал от М-оценок минимума максимальной асимптотической дисперсии параметров и получил [298] для функции потерь в крите рии (10.20) выражение (1 2) 2 при сout () =, (10.22) сout (1 2)сout при cout где cout – константа, зависящая от интенсивности промахов. Константу сout находят, численно решая уравнение [299] cout ) exp( 100 2 exp c 2 2.

= 2 d + out (10.23) 100 cout При 0.5 % 99 % зависимость сout() (рис. 10.2) удовлетворительно аппрокси мирует выражение cout = 1.211 exp(-/0.6642) + 1.656 exp(-/33.02). (10.24) Помимо (10.22), предложены [286, 291, 298, 364, 365] и другие способы задания, как полученные методом максимума правдоподобия, исходя из определенных требо ваний к оценкам, так и эвристические:

Lq – оценки: () = ||q, 1 q 2 (10.25) оценки Вэлша: () = d2 / 2[1-exp(-/d)2], d 0, (10.26) оценки Фэера: () = c2 [||/c - ln (1 + ||/c)], c 0, (10.27) оценки Демиденко: () = 2 / ( + d), d 0, (10.28) оценки Эндрюса:

(2) 1 {1 cos[(2 )1 / 2 ]}, | | (2 ) 1 / ( ) =, 0, (10.29) 1, | | (2) 1 / Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера 3. cout 2. 2. 1. 1. 0. 0 10 20 30 40 50 60 70 80 90, % Рис. 10.2. Зависимость параметра с out от интенсивности грубых промахов.

Точки – точные значения, линия – аппроксимация (10.24).

оценки Мешалкина:

() = -1 [1 - exp(-2/2)], 0, (10.30) оценки Рамсея:

() = -1 [1 - (1 + 1/2 ||)] exp(--1/2 ||), 0. (10.31) Последние три способа задания () предназначены для асимметричных функций распределения случайных погрешностей.

Все М-оценки обладают близкими свойствами, но считают [293], что предпочти тельны именно М-оценки Хьюбера, поскольку они сохраняют свои оптимальные свой ства при малом числе опытов и устойчивы и к несимметричным засорениям нормаль ного распределения погрешностей.

Критерий (10.20) – это гибрид критериальных функций, соответствующих зада чам МНК и метода наименьших модулей (МНМ). При интенсивности промахов константа сout, и функция М переходит в 1N.

U= (10.32) 2 k = При 100 % значение cout 0, и функция М становится эквивалентной критерию L (10.12). Подход Хьюбера «дал необходимое теоретическое обоснование методу наиме ньших квадратов и ряду других методов, приводящим к оценкам квазиправдоподобия.

При этом в вечном споре гауссовского и лапласовского распределений и соответствен но методов обработки по методу наименьших квадратов и по методу наименьших моду лей начинает просматриваться компромисс. Идеальным методом обработки признается тот, который сочетает в себе черты как одного, так и другого метода» [299].

Вычисляя М-оценки, измерения разбивают на две группы. Первая содержит экс периментальные точки, для которых |k| cout, вторая – остальные измерения, подвер женные, как считают, влиянию грубых выбросов.

Формула (10.20) справедлива, если Гаусс = 1. Поскольку М-оценки не инвари антны относительно преобразования масштаба, для обработки данных с Гаусс 1 пере ходят [298, 299] от (10.20) к другой минимизационной задаче:

N ~ [( k )] + a M = min.

M (, ) = (10.33) k = где коэффициент Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности ( )d, aM = (10.33) Nz () – функция плотности вероятности стандартного нормального распределения (0, Гаусс=1). В формуле (10.33) к неизвестным параметрам добавляется масштабный параметр. Величину 2 (дисперсию закона распределения погрешностей основных наблюдений Гаусс ) можно интерпретировать как выборочную остаточную дисперсию s0 для данных без промахов.

Задача (10.33) обобщается на многомерный случай, когда в каждой эксперимен тальной точке одновременно регистрируют свойств равновесной системы:

N ~ M (, ) = [( lk )] + a M = min, (10.34) l =1 k = ( ) d.

aM = (10.35) N z При использовании М-оценок Хьюбера пользователь должен варьировать про цент промахов, априорная информация о величине которого отсутствует. Тем самым процедура параметрической идентификации адаптируется к особенностям обрабаты ваемых данных.

Насколько эффективным будет анализ данных, зависит, не в последнюю очередь, от вычислительных свойств алгоритма М-оценивания и стратегии обработки реальных данных. Ниже в п. 10.3 описан разработанный нами вычислительный алгоритм CLINP 2.1 [366] робастной параметрической идентификации и обоснована общая стра тегия анализа первичных данных КФХА [367] (гл. 11).

10.3. Алгоритмы модифицированных методов Ньютона и Гаусса-Ньютона минимизации критериального функционала 10.3.1. Формулировка задачи Примем, что известны число сортов и стехиометрический состав комплексов, а подлежат расчету значения логарифмов констант устойчивости lg i, i =1, 2,..., p, и, возможно, факторы интенсивности li, l = 1, 2,..., ;

i = 1, 2, …,;

S, при общем числе искомых параметров z, находимом по формуле (9.127).

Точечные М-оценки параметров |lg i, li и – величины, соответствующие ~ ~ минимуму критериальной функции M. Критерий M можно минимизировать числен ными методами [1] непосредственно. Однако более эффективным считают [299] пе ~ реход от минимизации M к эквивалентной задаче метода взвешенных наименьших квадратов:

~ min M = + Q, (10.36) |,, где 1 N~ lk lk, = (10.37) 2 l =1 k = множители Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера | lk | сout, при ~ = lk lk (10.38) | lk | 2 lk сout, при 1 | | при lk c out ;

lk = (10.39) c out | lk | c out, при | lk | lk = w1 / 2 lk, lk = Alk - Alk, (10.40) lk wlk – статистические веса, назначенные в соответствии с моделью экспериментальных погрешностей, | | при lk cout ;

aM Q= (10.41) | lk | a M N cout / 2 при cout.

~ Чтобы избежать расчета производных M /, для минимизации (10.36) ис пользовали схему Хьюбера [298, 360]. Согласно этой схеме, располагая на m-й итера ции по оценкой (m), рассчитывают, не уточняя (m), неизвестные константы равнове сия и факторы интенсивности (вектор ), а затем по формуле (m) ( ) () 1N / ( m +1) 2 ( m) 2, kl = (10.42) ( m) a k =1 l = уточняют. Повторив расчет, заново вычисляют и т.д., до тех пор, пока изменение от итерации к итерации не становится меньше 1 %.

10.3.2. Расчет неизвестных факторов интенсивности Неизвестные параметры lg i и li рассчитывали по схеме Силлена (см. п. 9.2.2.4).

Как и в п. 9.6.1, выделили в отдельные подмножества реагенты с известными и неиз вестными факторами интенсивности li и представили матрицу А в виде (9.128). МНК оценки факторов интенсивности находят по формуле (9.130). Переход к М-оцениванию предполагает использование весов, и простое выражение типа (9.130) получить невоз можно. Чтобы воспользоваться М-оценками, преобразовали матрицы измеряемых свойств А, невязок и факторов интенсивности Е в такие векторы A#, #, Е#, что эле менты векторов yu соответствуют элементам матриц xlk при u = (l - 1) N+ k (l 1, 1 k N) для матриц А и ;

u = (l - 1) S + k (l 1, 1 k S) для матрицы Е. Так, элементам матрицы невязок элементы вектора # соответствуют следующим обра зом:

#1 = 11, #2 = 12,..., # = #+1 = 21, #+2 = 22,..., #2 =...

#1+(-1)N = 1, #2+(-1)N = 2,..., #N = N Затем сформировали весовую матрицу размером N с элементами lk = wlk lk, (10.43) Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности ~ а из нее – диагональную матрицу весовых множителей размером NN. Эле ~ менты последней = 1 / 2 при u = k + (l - 1)N. После выполненных преобразований uu lk равновесные концентрации реагентов размещаются в блочно-диагональной матрице С# размером SN, идентичные блоки которой – это матрицы С размером SN:

C 0... 0 0 C... C # = 0 0... 0............

0 0... C 1 1 В результате выполненных преобразований уравнения Фирордта можно перепи сать как A# = E# C#, (10.44) а формула (10.37) принимает вид 1 N # u ( u ) 2, # = (10.45) 2 u = или в векторном виде 1# [ ()]T # [ # ()], = (10.46) ~ ~ где # = T.

Теперь, чтобы найти неизвестные факторы интенсивности с использованием М-оценок Хьюбера, используем взвешенные переменные:

~ ~~ ~T ~ #T ~ ~ #~ A = A #, = #, C1 = C1, C1 = C1. (10.47) Окончательно формулу для расчета неизвестных факторов интенсивности получили в виде:

~ #~ ~T ~ # E1 = ( A E 2 C 2 ) C1 µ, (10.48) ~ ~T ~ ~T где µ = C1 (C1C1 ) 1.

~ 10.3.3. Производные взвешенных невязок по логарифмам констант устойчивости Для методов Ньютона и Гаусса-Ньютона необходимы первые и вторые производ ~ ные минимизируемой функции (10.36) по искомым параметрам M / i = Ф/ i, ~ 2 M / ij = 2Ф/ ij. Учитывая, что факторы интенсивности вычисляются по опи санной выше схеме, а для итеративного уточнения производные не требуются, необ ходимы лишь производные по параметрам lg i. Для их расчета, в свою очередь, тре ~ ~ буются производные u / lg i и 2 u / lg i lg j. Начнем с частного случая, когда w в минимизируемом функционале M (10.36) веса lk = 1, а расчету подлежат производ ные lk / lg i, 2lk / lg ilg j. Введя обозначения Di = / lg i, и Dij = 2 / ( lg i lg j ), из формулы (9.137) получаем:

{ } ( Di C1T )µ T + E ( Di C ) (1 µC1 ), Di = Di A = (10.49) ln Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера { } ( Di A)( D j C1T ) + ( D j A)( Di C 1 T ) + ( Dij C 1 T ) µ T + 2 Dij = (ln 10) (10.50) { } 1 + E ( Dij C ) + ( Di E 1 )( D j C 1 ) + ( D j E 1 )( Di C 1 ) (1 µC 1 ).

(ln 10) В общем случае неединичных весов (lk 1) применяем взвешенные варианты ~~~ ~ ~ переменных A,, C1, C1T, µ и приходим к следующим выражениям:

{ }, ~ ~T ~ ~ ~~ ~ ~ ( Di C1 )µ T + E ( Di C )(1 µC1 ) Di = Di A = (10.51) ln {(Di A)(D j C1T ) + (D j A)(Di C1T ) + (Dij2C1T )} µ T = 2~ ~ ~ ~ ~ ~ ~ 2~ ~ Dij = Dij A = (ln 10) (10.52) { } 1 2~ ~ ~ ~~ =+ E ( Dij C ) + ( Di E1 )( D j C1 ) + ( D j E1 )( Di C1 ) (1 µC1 ).

(ln 10) 10.3.4. Итеративная вычислительная схема методов Ньютона и Гаусса-Ньютона Для расчета lg применяли модифицированные методы Ньютона и Гаусса-Нью тона.

~ Пусть в окрестности точки lg c функцию M можно аппроксимировать квадра тичной формой ~ ~ M (lg ) M (lg c ) + g (lg c ) T (lg lg c ) +. (10.53) + (lg lg c ) T H (lg c ) (lg lg c ), где градиент ~ g = M = = J (lg ) T # # (lg ), (10.54) J(lg ) – матрица Якоби c элементами ~ Jui(lg ) = u / lg i, u = 1, 2,..., N, i = 1, 2,... p, (10.55) матрица Гессе ~ Н = 2 M = 2 = J(lg )T # J(lg ) + F, (10.56) (J (lg ) T ) # # (lg ), F= (10.57) ~ [J( lg ) ]uij = 2 u / lg i lg j = [Jui(lg )] / lgj. (10.58) В методе Ньютона матрицу Гессе вычисляли точно, а в методе Гаусса-Ньютона элементы матрицы F приравнивали нулю.

Из формулы (10.53) следует (см. п. 9.2.2.4) схема итеративного уточнения lg c:

lg + := lg c + s, (10.59) + где s – направление движения от текущей оценки lg c к уточненной lg, s = -H-1 g = -J(lg c)T # #( lg c). (10.60) При использовании метода Ньютона s = { J ( lg c ) T # J ( lg c ) + ( J ( lg c )) T # # ( lg c )} (10.61) J ( lg c ) T # # ( lg c ), Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности а в методе Гаусса-Ньютона расчет s упрощен:

( )1 J (lg c ) T # # (lg c ).

s = J (lg c ) T # J (lg c ) (10.62) Из-за приближенного расчета матрицы Гессе метод Гаусса-Ньютона не приводит в точку минимума, если (а) из-за неудачного выбора химических форм, учтенных в мо ~ дели, и/или плохого выбора начальных приближений параметров велики невязки u и ~ (б) функция M (lg ) существенно нелинейна (велики вторые производные J (lg ) ) [328]. В то же время, алгоритм метода Гаусса-Ньютона экономно потребляет ресурсы ЭВМ и сходится к минимуму в подавляющем большинстве задач.

10.3.5. Обеспечение глобальной сходимости методов Ньютона и Гаусса-Ньютона Как указано в п. 9.2.2.4, для обеспечения глобальной сходимости и высокой ско рости расчетов классическую вычислительную схему метода Ньютона необходимо мо дифицировать.

Проблеме глобальной сходимости методов ньютоновского типа посвящены мно гочисленные монографические и оригинальные работы (см., например, [323, 327, 328, 369–373]). В нашем алгоритме CLINP 2.1 воплощены результаты исследований пос леднего времени.

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

Для коррекции неположительно определенной матрицы Гессе использовали ме тод, основанный на ее спектральном разложении [323, 328, 374]. Исходную матрицу Н представляли в виде произведения трех матриц H = U UT, (10.63) где К – диагональная, U – ортогональная матрицы, составленные из собственных чисел ~ и собственных векторов матрицы Н. Модифицированную матрицу H получали, заме няя отрицательные собственные значения i их абсолютными величинами:

~ ~ H = U K UT, (10.64) ~ где K – диагональная матрица с элементами ~ ii = ii при ii 10 16, ~ = ~ = | | при ii 10 16, ii ii (10.65) ii ~ при 10 16 ii 10 16.

= ii При такой модификации матрицы Гессе вектор поправок s, вычисленный с использова ~ нием матрицы H, имеет ту же норму, что и вычисленный обычным методом.

В ходе итераций алгоритм может попасть в седловую точку lg c. В седловой точке градиент g(lg c) равен нулю, вследствие чего и вектор поправок s будет нуле вым. Чтобы обойти седловую точку, надо выполнить шаг в направлении «отрицатель ной кривизны»2 [323, 328, 372]. Таким направлением будет любой вектор d, удовлетво ряющий неравенству В методе Гаусса-Ньютона матрица Гессе всегда определена неотрицательно, а при матрице Якоби полного ранга – положительно [286].

Из-за упрощенной оценки матрицы Гессе в методе Гаусса-Ньютона найти в этом методе направление отри цательной кривизны не удается. Как результат, в окрестности седловой точки итерации сходятся именно к ней.

Глава 10. Модель химических реакций: оценивание параметров при нарушении стандартных предпосылок МНК. М-оценки Хьюбера dT H d 0. (10.66) В алгоритме CLINP 2.1 в качестве направления d выбран собственный вектор, соответствующий наибольшему по абсолютной величине отрицательному собствен ному значению матрицы Н.

Таким образом, для движения из точки lg c предлагаются два направления: нью тоновский (квазиньютоновский) шаг s и направление отрицательной кривизны d. Эти векторы называют убывающей парой («descent pair») [328]. На ее основе строили шаг, обеспечивающий глобальную сходимость алгоритма и квадратичную скорость ло кальной сходимости.

К текущим значениям lg c на каждой итерации вводили поправки :

lg +() = lg с + (), (10.67) где 0 1 – длина шага, вектор – линейная комбинация векторов s и d (в методе Гаусса-Ньютона d = 0). Следуя [373], вектор () находили как () = s + d. (10.68) При движении вдоль вычисленного таким образом направления критериальная функция убывает, но при большой длине шага можно «проскочить» минимум и по лучить такие lg +, при которых критериальная функция возрастет. Одним из надежных способов обезопасить себя от подобного поведения алгоритма является метод линей ного поиска длины шага [327, 328]. При линейном поиске обеспечено движение к ми нимуму из начальных точек даже вне области сходимости метода Ньютона.

Известно [369, 370], что для сходимости итераций к минимуму необходимо вы полнение неравенств1:

Ф(lg с+)Ф(lg c)+Ф(lg c)T, (0, 1), (10.69) или Ф(lg +)T Ф(lg c)T, (, 1), (10.70) где и – некоторые константы. Достаточным условием сходимости алгоритма явля ется соблюдение обоих неравенств (10.69) и (10.70) [375]. Выполнение условия Армио (10.69) гарантирует достаточно ощутимое убывание Ф по сравнению с длиной шагов, а условия (10.70) – достаточно большую длину шагов.

Для того чтобы избежать чрезмерно малых шагов, величину ограничили снизу значением мин;

чтобы не допустить выхода алгоритма из области сходимости, макси мальный шаг ограничили величиной hмакс. При |||| (hмакс)2 максимальную длину шага находили как макс = hмакс / ||||1/2. (10.71) Таким образом, длину шага выбирали из интервала [мин;

макс] по следующей схеме [328]:

ШАГ 1. Положить = 1.

ШАГ 2. Найти |||| и назначить при || || (hмакс ) 2, макс = (10.72) hмакс / || || 1 / 2 при || || ( hмакс ) 2.

ШАГ 3. Оценить Ф(макс) = Ф(lg с + макс). (10.73) Если Ф(lg с + макс) удовлетворяет условию (10.69), то принять макс как длину шага. В противном случае перейти к шагу 4.

Соблюдение условия Ф(lg + = lg с + ) Ф(lg с) не всегда обеспечивает сходимость итераций к точке минимума.

Часть III. Вычислительные аспекты количественного физико-химического анализа. Построение моделей и оценка их достоверности ШАГ 4. С использованием Ф(макс), Ф(=0) = Ф(lg с) и Ф/(=0) = Ф(lg с)T, (10.74) аппроксимировать Ф() квадратичной функцией mq() = [Ф(lg с + макс) - Ф(lg с) - Ф(lg с)T ]2 +Ф(xс)T + + Ф(lg с) = [Ф(=макс) - Ф(=0) - Ф(=0)]2 + Ф(=0) + Ф(=0). (10.75) ШАГ 5. Рассчитать значение (*), минимизирующее mq():

(*) = - Ф(=0) / (2{Ф(=макс) - Ф(=0) - Ф(0)}). (10.76) (1) ШАГ 6. Выбрать длину шага как (1) = max[мин;



Pages:     | 1 |   ...   | 3 | 4 || 6 | 7 |   ...   | 10 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.