Полак__Применение_вычислительной_математики_в_химической_и_физической_кинетике (972294), страница 19
Текст из файла (страница 19)
Из формул (4), (5) следует, что максимуму вероятности соответствует максимальное значение функции правдоподобия. В прикладных задачах обычно максимизируют 1п Ь (с, 6), что более удобно. Часто принимается 137), что получаемые опытным путем данные распределены по нормальному закону. В таком случае будем иметь где с„(0) — концентрация реагента, вычисленная для условий и-го опыта при заданных значениях параметров 0; а,~ — дисперсия, характеризующая рассеяние экспериментальных данных. Учитывая выражение (5), получаем следующее выражение для функции правдоподобия: 1п Ь(с, О) = сопз[, — — ~ с„'[с — с„(0)['. (7) Из (7) следует, что максимуму правдоподобия соответствует минимум взвешенной суммы квадратов отклонений вычисленных значений концентраций от опытных, т.
е. принцип Фишера сводится к известному методу наименьших квадратов. В качестве весов служат обратные значения дисперсий. Так как почти всегда дисперсии неизвестны, их приходится заменять выборочными значениями г~. В этом случае плотность распределения опытных данных будет характеризоваться законом Стьюдента [33). Функция правдоподобия представится в виде и 1пЛ(е 0) сопзЬ вЂ” — ~~~~~ (~„+1)1п 1+ [с„— с,„(0) [~, и=1 ~ и~и (8) где ~ — число степеней свободы выборочной дисперсии г'. В практической работе, однако, общепринято пользоваться методом наименьших квадратов и в тех случаях, когда отсутствуют сведения не только о теоретических (генеральных) дисперсиях аз, но и о выборочных з'. Насколько оправдан такой подход, неизвестно.
Между тем степень приближения можно оценить, воспользовавшись расчетом вероятностных моделей некоторых простых случаев (например реакции первого порядка). Несмотря на определенный интерес данного вопроса, он до сих кор никем не рассматривался. Следует заметить, однако, что если справедливо предположение о нормальном законе распределения опытных концентраций, то оценки кинетических параметров, получаемые максимизацией выражения (7), будут несмещенными, достаточными и асимптотически нормальными и в тех случаях, когда вместо оз в (7) подставляются значения г'. Что касается дисперсии определяемых параметров, то, исходя из общей теории максимального правдоподобия [33), можно утверждать, что оценка неизвестных параметров 0 по формуле (8) будет более точна, чем по (7).
В кинетических исследованиях диапазон превращений реагирующих веществ обычно стараются поддерживать довольно большим [31[. При атом концентрации отдельных компонентов могут изменяться в ходе реакции в несколько раз. Трудно ожидать, чтобы в этих случаях абсолютная величина рассеяния экспериментальных данных в каждой точке была одинакова. Опытные данные как правило, свидетельствуют о том, что остается примерно одинаковой относительная величина ошибки, в то время как абсолютное ее значение зависит от величины концентрации.
На рис.28 приведен график, построенный по данным, полученным Фейгиным при изучении пиролиза метана в трубчатом реакторе [20, 381, который отражает подобную зависимость. Как бснй ~о' з5 Рнс. 28. Эазпспмость среднекзадратнческой ошцбкп определения концентрацпп метана г(СН4) ог велнчпны его концентрапнн [20, 38! 5 и И [Сиг), абылл:% видно из рис. 28, величина среднеквадратической ошибки увеличивается с ростом концентрации реагирующего вещества. Отношение з (с)/с (где с — концентрация, г (с) — ее ошибка) принято называть коэффициентом вариации [37). Если коэффициент вариации постоянен, это указывает на наличие логарифмически нормального распределения [32, 37, 39]. При этом максимуму правдоподобия будет соответствовать минимум суммы квадратов отклонений логарифмов концентраций. Однако при изучении кинетики пиролиза метана применялся обычный метод наименьших квадратов [20).
В атой же работе приводятся также результаты других авторов, применявших при обработке данных по пиролизу СНа логарифмически нормальное распределение. В этом случае удалось получить кинетические параметры, более близко отвечающие теории. Учитывая зависимость ь (СН,) от [СНа), приведенную на рис. 28, такой вывод не представляется неожиданным. Сказанное здесь подтверждается также работой [24), где для нахождения наиболее вероятных значений констант применялся логарифмический вариант метода наименьших квадратов, так как было показано постоянство коэффициента вариации.
В литературе по теории вероятностей [401 подчеркивается, что использование «нормальной техники» в случае логарифмически нормального аакона распределения может привести к существенным ошибкам, особенно при обработке небольшого числа экспериментов. В статье [41) приводится перечень некоторых критериев, соответствующих различным распределениям, которые могут встретиться в экспериментальной работе.
Другие примеры применения принципа максимального правдоподобия для определения констант скоростей химических реакций приведены в работах [42, 43). й. Локалъмьае мепаоды таомсма монсэманпа еыороетпей После того как выбран критерий оптимальности кинетических параметров, зависящий, как подчеркивалось выше, от вида распределения опытных значений концентраций, можно перейти к определению таких величин параметров, которые обращают выб: ранный критерий в минимум (или максимум). В общем случае может быть несколько экстремумов критерия.
Задача будет заключаться в отыскании координат экстремума, наибольшего по абсолютной величине. Для решения таких задач разработаны различные локальные и нелокальные методы поиска (см. обзорные статьи [7, 44, 45[ и книги [46,47[). Однако пока еще не все эти методы нашли применение для обработки кинетических данных. Поэтому здесь будут рассмотрены лишь наиболее распространенные методы поиска констант — вначале локальные, затем нелокальные.
Для удобства изложения мы будем предполагать, что опытные данные распределены по нормальному закону с одинаковой дисперсией. В этом случае принцип максимума правдоподобия сводится к методу наименьших квадратов. Наиболее вероятными будут такие значения параметров, которые минимизируют сумму квадратов отклонений вычисленных величин концентраций от их опытных аналогов * Я(О) =,'Я [с„— с„(0))', (9) где Я (О) — сумма квадратов отклонении; с„— концентрация компонента в и-м опыте (и = 1, 2, ..., Дт); с„(0) — вычисленное значение концентрации для условий и-го опыта, зависящее от вектора кинетических параметров О = (0~), где [ = — 1, 2, ..., р. Поскольку схема реакции и опытные данные являются заданными, с„(0) и Я (О) будут функциями только кинетических параметров О.
Рассмотрим теперь методы определения локального минимума функции Я (0), начав с градиентных методов. а) Метод крутого спуска Градиентные методы в отличие от простых методов, описанных выше, основаны на использовании сведений о наклоне поверхности критерия в данной точке. Движение к минимуму суммы квадратов отклонений производится по направлению наибольшей а Кслп в эксперпменте известны концентрации нескольких веществ, то суммпрованне в выражении (9) необходимо выполнить для всех этих веществ. крутизны поверхности Я (О), которое противоположно направлению градиента: р 8гад Ю(0) =,Я „( ООо где 00; — единичные векторы в направлении координатных осей. Производные дЯ (О)/дО; в формуле (10) обычно находятся численно, для чего каждый из параметров О; (1 = 1, 2, ..., р) меняется поочередно на некоторую небольшую величину (например ЛО; = == (0,05 — 0,1)0;).
Производные заменяются величинами ЬЯ (0)/ЛО;. Как показано Лапидусом и сотр. [48], такой подход к определению проиаводных неэффективен, поскольку параметры варьируются лишь в одну сторону. Авторы [48] предлагают вычислять производные с помощью так называемых факторных планов (см., например, [49]). Более точное определение производных дЯ (О)/дО~ позволяет сократить общее число необходимых вычислений [48] Можно найти эти производные иначе: = — 2 Я~ [си — сн(0)]лги п=а (1.. 1,2,...,р), (11) где дс„(э) (12) дэ; Если кинетика исследуемых реакций описывается системой обыкновенных дифференциальных уравнений, то эту систему можно дополнить дифференциальными уравнениями для производных ' — „"'," =-/го ('.(0), О] (и=1,2,...,Х; 1=1,2,..., р) (13) с начальными условиями х; = О [50].
Правые части уравнений (13) находят согласно правилам, указанным в работе [50] и подробно описанным на стр. 66. Интегрируя (13) совместно с уравнениями кинетики, мы получаем непосредственно значения производных де„(0)/дО;. Подставляя их в вырансение (11), а (11) в (19), находим направление и величину градиента функции Я (О).
Такой способ вычисления производных применялся в ряде работ по изучению кинетики сложных химических реакций [10, 51 — 53]. Применение факторных планов для этих целей описано в работе Бокса [54] иа примере двух последовательных мономолекулярных реакций ьз А — В-' С. В этой работе, а также в книге [491 дается описание последовательности вычислений, необходимых для движения по линии крутого спуска. Последнее производится до тех пор, пока Я (О) снижается.