Sensitivity (972306), страница 6
Текст из файла (страница 6)
ясно, что константу k2 невозможно определить из кинетической кривой вещества A, так как скорость расхода A определяется только константой k1. Если же модель реакции насчитывает десятки элементарных стадий, то подобные факты далеко не очевидны.
В данной работе была изучена возможность применения коэффициентов чувствительности для отбора экспериментальных данных. Исследование проводилось на простом модельном примере последовательной реакции с двумя наборами констант скорости: а) k1k2; б) k1, k2. В обоих случаях использовались начальные условия A0, B0C0. Задача заключалась в выявлении корреляции между величиной коэффициентов чувствительности и точностью определения соответствующих констант скорости.
Экспериментальные данные имитировались путем добавления случайных ошибок с нормальным распределением и заданной дисперсией к точным решениям прямой кинетической задачи с известными константами скорости. Для генерации данных была составлена программа, использующая алгоритм Бокса-Мюллера-Белла [18] для преобразования случайных чисел с равномерным распределением на интервале в случайные числа с нормальным распределением, средним значением 0 и дисперсией 1. Случайные числа с равномерным распределением получались от стандартной подпрограммы RANDOM.
Для каждой кинетической кривой было построено по 5 наборов данных, моделировавших серии независимых измерений. Во всех случаях уровень «экспериментальных» ошибок составлял 10% от среднего значения концентрации на рассматриваемом участке кинетической кривой. Для решения обратной задачи выбирали тем или иным способом 8 точек на одной из кривых и определяли одну или две константы скорости.
Для первой задачи (k1k2) исследовалась надежность определения k1 и k2 из разных участков кинетической кривой промежуточного продукта B. С точки зрения анализа чувствительности можно выделить две области на кривой: область малых времен (t1.6) и область больших времен (t1.6). В первой области коэффициенты чувствительности B к константам k1 и k2 приблизительно одинаковы. Во второй области чувствительность к k1 примерно в 4 раза ниже по сравнению с чувствительностью к k2 (см. рис. 1). В соответствии с этим были сформированы две группы данных — 5 независимых серий по 8 точек с шагом 0.2 в интервале 0.2t1.6 и 5 таких же серий в интервале 1.6t3.0 с. Решая обратную кинетическую задачу, мы получили по 5 оценок каждой константы для двух областей кинетической кривой. Далее по формулам элементарной статистики были вычислены средние значения и дисперсии найденных оценок. Результаты приведены в таблице 10. Первые 5 строк содержат оценки констант для 5 серий модельных данных. В скобках указаны оценки стандартных отклонений, полученные из ковариационной матрицы, которая вычислялась в приближении Гаусса-Ньютона, как описано в следующем разделе. Последняя строка содержит средние значения и стандартные отклонения, полученные по 5 оценкам. Хотя средние значения констант скорости во всех случаях оказались достаточно близкими к истинным значениям, дисперсии обнаруживают явную корреляцию с чувствительностями. Как следует из таблицы, в области малых времен обе константы определяются приблизительно с одинаковой степенью надежности, тогда как во второй области оценка k2 существенно надежнее, чем k1.
Для второй задачи (k1, k2) константы скорости определялись из кинетической кривой конечного продукта C. В этом случае мы имеем дело с различием нормированных коэффициентов чувствительности на 2 порядка практически вдоль всей кинетической кривой. Так же, как и в первой задаче, были обработаны 5 независимых серий данных по 8 точек в каждой. Точки были выбраны в интервале 0.6t2.0 с. Полученные результаты представлены в таблице 11. В то время как k1 определяется с той же степенью надежности, что и в первой задаче, получить оценку k2, по существу, не удалось (стандартное отклонение в два раза превышает оцениваемую величину). Следует также отметить, что для одной серии данных программа вообще не смогла довести до конца минимизацию отклонений, так как значение k2 оказалось слишком большим, и концентрация C стала практически независимой от дальнейшей вариации этой константы.
Таким образом, рассмотренные модельные задачи иллюстрируют ситуации, когда для получения оптимальных оценок констант скорости важен выбор подходящей кинетической кривой, а также участка на этой кривой. Локальный анализ чувствительности во всех случаях оказывается руководством к правильному выбору.
Таблица 10. Результаты определения констант скорости k1 и k2 по двум участкам кинетической кривой B.
0.2t1.6 | 1.6t3.0 | ||
k1 | k2 | k1 | k2 |
| | | |
| | | |
| | | |
| | | |
| | | |
| | | |
Таблица 11. Нормированные коэффициенты чувствительности и результаты определения констант скорости по кинетической кривой [C].
III.2. Минимизация отклонений и коэффициенты чувствительности
Известные методы поиска минимума можно разделить на три группы. В первую группу попадают методы, основанные на вычислении лишь значений минимизируемой функции. Вторая группа образована методами, требующими вычисления первых производных (градиента) минимизируемой функции. Наконец, в третью группу попадают методы ньютоновского типа, для применения которых необходимо уметь вычислять матрицу вторых производных (матрицу Гессе) рассматриваемой функции. Вообще говоря, методы ньютоновского типа наиболее эффективны (требуют минимального числа итераций), но их применению часто мешает сложность получения вторых производных. Градиентные методы также весьма эффективны и широко используются, поскольку первые производные обычно проще получить, чем вторые.
Пусть ищется минимум функционала
относительно варьируемых констант скорости k1,,kM, где величины i представляют собой разности рассчитанных и заданных концентраций веществ в различные моменты времени:
(индекс i нумерует экспериментальные точки, которые, вообще говоря, отвечают разным веществам и разным моментам времени, т. е. i заменяет комбинации индексов n и m; L – общее количество экспериментальных данных).
Дифференцируя (24) по константам скорости, получим элементы вектора градиента:
Очевидно, входящие в правую часть уравнения (25) частные производные есть не что иное, как коэффициенты чувствительности веществ n по отношению к константам скорости kj, вычисленные для моментов времени tm. Аналогично можно получить выражения элементов матрицы Гессе функции F через коэффициенты чувствительности 2‑го порядка.
Таким образом, умение вычислять коэффициенты чувствительности 1‑го или 2‑го порядка позволяет реализовать эффективную процедуру решения обратной кинетической задачи. Однако, заранее трудно сказать, можно ли при этом достичь серьезного выигрыша в общем объеме вычислительной работы. При переходе к более быстро сходящемуся методу градиентного или ньютоновского типа уменьшится число итераций, но возрастут затраты на одну итерацию, поскольку вместо одной системы дифференциальных уравнений придется решать N систему того же размера. При этом, как было показано в п. II.3, наблюдается дополнительный рост затрат пропорционально количеству моментов времени tm, для которых имеются данные.
В программе KINET обратная кинетическая задача решается с помощью метода Пауэлла [31‑33], рассчитанного специально на минимизацию суммы квадратов и не требующего вычисления производных. Мы попытались оценить, возможен ли выигрыш в случае перехода к более быстро сходящемуся методу с применением коэффициентов чувствительности.
В таблице 12 представлены величины, характеризующие объем вычислений в процессе решения двух простых обратных задач. В первом случае по 4‑м экспериментальным точкам оценивалось значение одной константы скорости, что потребовало 22 раза решить прямую задачу. Во втором случае 2 константы скорости определялись по 8 точкам на кинетической кривой одного вещества, что потребовало 67 раз решить прямую задачу. В той же таблице приведены затраты на однократное вычисление набора коэффициентов чувствительности. Такой объем работы потребуется на каждой итерации, если будет использована процедура минимизации с вычислением производных. Предположим, что мы реализуем метод наименьших квадратов Гаусса-Ньютона [3]. Это один из наиболее быстро сходящихся методов минимизации суммы квадратов отклонений, обладающий свойствами ньютоновского метода, но требующий вычисления лишь первых производных. Как правило, при наличии сходимости метод Гаусса-Ньютона дает результат за 4‑6 итераций. Сопоставляя это с данными таблицы 12, можно сделать вывод, что при переходе от метода Пауэлла к методу Гаусса-Ньютона объем вычислений в лучшем случае останется на прежнем уровне, и никакого выигрыша не ожидается. Конечно, полученная нами оценка касается лишь рассмотренных простых задач. Не исключено, что при переходе к более сложным задачам соотношение затрат может измениться в пользу метода с вычислением производных.
Таблица 12. Оценка вычислительных затрат при решении обратной кинетической задачи.
Последовательная реакция A B C, определение k1 | |||
Вычисление минимизируемой функции – 22 раза | |||
Решение прямой задачи до t 0.8 сек: | |||
Число шагов | Вычисление пр. частей ДУ | Вычисление якобиана | |
Однократно | 50 | 89 | 6 |
Всего | 1124 | 2206 | 123 |
Вычисление матрицы ковариаций (расчет функций Грина) | |||
Число шагов | Вычисление пр. частей ДУ | Вычисление якобиана | |
303 | 677 | 790 | |
Последовательная реакция A B C, определение k1 и k2 | |||
Вычисление минимизируемой функции – 67 раз | |||
Решение прямой задачи до t 3.0 сек: | |||
Число шагов | Вычисление пр. частей ДУ | Вычисление якобиана | |
Однократно | 64 | 143 | 6 |
Всего | 4200 | 9269 | 407 |
Вычисление матрицы ковариаций (расчет функций Грина) | |||
Число шагов | Вычисление пр. частей ДУ | Вычисление якобиана | |
655 | 1428 | 1647 |
III.3. Получение матрицы ковариаций для констант скорости.
При решении обратной кинетической задачи важно не только оценить значение неизвестной константы скорости, но и охарактеризовать надежность полученного результата величиной дисперсии, стандартной ошибки или доверительного интервала. В задачах нелинейного МНК, к которым относится и обратная кинетическая задача, вычисление указанных характеристик обычно проводят по тем же формулам, что и в линейном МНК [1, 3]. Такой подход является приближенным и основан на уже упоминавшемся выше методе Гаусса-Ньютона. Условием минимума функционала (24) является равенство нулю вектора градиента, элементы которого даны выражениями (25). Это условие может быть представлено матричным уравнением