Оценивание параметров обыкновенных дифференциальных уравнений 653
Примеры оценивания, использующего нелинейные (по пара метрам) аналитические решения детерминированной модели про цесса, соответствующей выражению (9.1.7), уже приводились ранее в гл. 6. Поэтому ниже метод наименьших квадратов будет приме
нен к модели процесса, нелинейной |
по зависимым переменным |
[аналогичной уравнению |
|
|
(9.1.10)] и, |
следовательно, |
План |
эксперимента |
требующей |
численного ре |
шения. |
|
|
(гл.8) |
|
|
|
Если для решения диф ференциальных уравнений используются разностные уравнения, то весьма целе сообразно сочетать стан дартную программу нели нейного оценивания с под ходящей стандартной вы числительной программой. Многошаговый метод Адамса — Мултона и одношаговый метод Рунге — Кутта хорошо известны как стан дартные методы решения задач с начальными усло виями [5, 6]. Большинство вычислительных центров
Ф и г. 9.2.2. Поток информа ции при оценивании парамет ров модели с начальными зна чениями с помощью числовых разностных схем.
Сбор данЧЫХ:
вход X, выхоі1Y
Введение данных и модели |
в машинную |
программу |
Выбор предполагаемых значений параметров
Применение разностной схембЛ [длярешения дифференциально гоуравнения для каждого зле -
мента »х
Нелинейное оценивание |
параметров |
и соответствцющих |
статистик |
(гл. В) |
|
Анализ (гл.5,6 и 7)
имеет вычислительные коды для получения решений. С помощью аналоговых и гибридного типа вычислительных машин или цифро вых модуляторов аналоговых вычислительных машин (типа МіМіС) могут легко выполнять расчеты по стандартным кодам даже лица, плохо знакомые с программированием. Прежде всего необходимо за дать начальные оценки коэффициентов модели, которые затем уточ няются с помощью некоторой итерационной схемы, представленной в общих чертах диаграммой потока информации на фиг. 9.2.2. Конечно, в численных схемах и при округлении могут встретиться трудности. В числовых разностных схемах Рунге — Кутта и дру гих может возникнуть числовая неустойчивость, если нарушаются
некоторые критерии устойчивости [7]. В некоторых моделях боль шие члены приблизительно равной величины вычитаются друг из друга, что приводит в результате к потере значащих цифр.
Пример 9.2.1. Оценивание кинетических констант скорости по дискретным наблюдениям
В работе [8] была применена процедура, намеченная в общих чертах на фиг. 9.2.2, для оценивания кинетических констант про цесса, который осуществляется с участием некоторых металлоорганических соединений по следующей схеме:
ki
1 + 3 ^ 2 + 2,
2 + 4 7^3 + 3, hs
3 + 5 ^ 4 + 4, ki
1 + 4 ^ 2 + 3,
2 + 5 ^ 3 + 4,
ke
1 + 5 ^ 2 + 4,
где kj — константы скоростей прямых реакций, а цифры обозна чают химические вещества. Сначала предполагалось, а затем было проверено, что константы скоростей не зависят от состава смеси. Шесть констант равновесия для реакций (КІ) вычислялись отдель но, так что требовалось оценить лишь половину из двенадцати констант скоростей (константу скорости обратной реакции можно рассчитать по константе скорости прямой реакции).
Из механизма реакции следует, что реагируют пять веществ 1—5, которые подчиняются двум уравнениям полного баланса вещества, опирающимся на первоначальный состав при t = О, вследствие чего для построения модели требуется лишь три неза висимых дифференциальных уравнения. Оставшиеся компоненты затем можно было бы получить из полного баланса вещества. Эти ми тремя уравнениями были следующие:
^ _ 2 * . ( • , * - • £ . ) - Ъ |
+*« |
- |
Оценивание |
параметров обыкновенных дифференциальных |
уравнений |
655 |
+ Ä 4 ( с і С 4 — + & 5 ( с 2 С 5 — - ^ - ) »
Начальными условиями служили известные детерминированные
значения: Cj (0) = |
с1 0 ; с2 (0) = |
0; с3 (0) = 0; с4 (0) = |
0 и с5 (0) = |
= с5 0 . Схема интегрирования |
была взята из работы [9], где пред |
лагался ряд быстро |
сходящихся и точных формул типа |
прогноз — |
коррекция. Конкретный метод сводился к продолжению решения
уравнения old dt — С = i (t, С) |
от |
значений, уже |
найденных |
в предыдущие моменты времени t0, |
ti, |
. . ., tn. Пусть Cn |
— вектор |
значений состава в момент tn. Тогда предлагалась следующая про
цедура {h = |
tn+i — |
tn): |
1. Предсказывается (оценивается) следующее значение векто |
ра зависимой |
переменной: |
|
Си+і = |
Сп_з - j - -g- (2Cn — C„_i - j - 2Cn _2 ). |
2.Вычисляется производная
3.Уточняется оцененное значение
C n + i = Ï7 ~Ь Ï7 ^ " - i — Ï7 |
~Ь ^ ( ï 7 |
~b f f ^ n ) • |
4. Вычисляется производная
Cf'i+i — i (^n+ii Cn +i).
5. Проверяется ошибка усечения с использованием значения (27/503) (С£+і — С „ + і ) . Как указали Болл и Гроенвеге, уравнения Фельберга [9] требуют меньше машинного времени, чем фор мулы интегрирования Хэмминга [10].
Для оценивания методом наименьших квадратов использовал ся метод Маркуардта, описанный в гл. 6, в частности выражение (6.2.20). Необходимые частные производные получались численно путем интегрирования уравнений для скоростей реакций со слегка возмущенными константами скоростей и вычисления соответствую-
656 Глава 9
щих |
разностей |
|
д £ ч |
^Cjjjkj, |
. . . , ki + Mti, |
kL) — Cj](ki |
дкі |
|
|
Акі |
Следовательно, для получения необходимого числа значений про изводной одна итерация Маркуардта требовала по крайней мере семикратного интегрирования уравнений для скоростей реакций.
В одном из приложений описанной выше программы измеря лось содержание каждого из пяти компонентов (в г-моль/л) в следу-
|
|
1,00 |
|
|
|
|
\ |
|
|
1 |
|
|
0,80 |
- |
|
|
|
|
|
|
|
|
« |
0,60 L |
|
|
Ci |
|
A |
|
|
|
I |
• |
|
|
|
|
|
|
|
|
Y |
|
|
С3 |
|
|
|
|
|
3 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
S 0,20 ТА |
• |
|
|
|
|
—tr- |
L |
|
j a |
0,00 Р Ѵ * ~ |
' |
І |
|
1 |
і. |
|
|
С г |
A— |
f—1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
3 |
1000 |
|
2000 |
|
|
3000 |
|
|
|
Время |
реакции, ми» |
|
|
|
|
Ф и г. П . 9 . 2 . 1 . |
Рассчитанные |
кривые |
откликов и экспериментальные дан |
|
|
|
|
|
ные |
[8]. |
|
|
|
|
9 |
|
|
Экспериментальные |
точки: |
|
|
(трет-С4 Нв О)4 Ті;/Д (mpem-C4H.O)3[(CH.,)2N] Ті; |
О (mpem-C4 H„0)2 [(CHa )jNl2 Ti; |
|
|
Д (mpem-C«H,0)[(CH3)2N],Ti; |
• C(CH3 )2 N]4 Ti. |
|
|
|
Кинетические |
константы: |
|
|
fei = |
3,74 - Ю - 3 , fc2"= 8,33-10-», ft, |
= 5,33-10-«, ft4 |
= 2,47-10-*, |
ft, = 1,58• 10-«, ft, = |
|
|
|
|
= |
2 , 5 8 - Ю - 5 . |
|
|
|
|
|
|
|
Константы |
равновесия: |
|
|
|
|
Кі |
= 7,40 - Ю - 2 , |
Кг |
= 2 , 9 7 - Ю - 1 , |
К, |
= |
9,50 10-«. |
ющие моменты времени (в минутах): 0, 10, 40, 70, 202, 490, 1190, 1453, 2410, 2795 и 3765. Таким образом, для оценивания шести констант скоростей использовалось полное число 5-10 = 50 изме ренных значений состава. Исходные оценки констант скоростей получались подсчетом вручную тангенса угла наклона в несколь ких точках на кривых зависимости концентрации от времени. Значения функции ф на каждой итерации оказались следующими:
Итерация |
^ |
IBM 7040 |
Машинное время, с |
0 |
|
0 |
1 |
0,0278 |
100 |
2 |
0,0254 |
186 |
3 |
0,0215 |
273 |
4 |
0,0206 |
360 |
Оценивание |
параметров |
обыкновенных |
дифференциальных |
уравнений |
657 |
Константы скоростей изменялись
Начальные оценки |
Конечные оценки |
, значений h |
значений ft |
kt: |
0,053 |
0,0418 |
k2: |
0,025 |
0,0272 |
k3: |
0,0048 |
0,0218 |
следующим образом:
Начальные оценки |
Конечные оценки |
значений h |
значений ft |
А4 : |
0,0085 |
0,0208 |
кь: |
0,0069 |
0,000532 |
кв: |
0,0141 |
0,0108 |
Оценки быстро сходились при очень малом числе итераций, если предполагаемая кинетическая модель хорошо описывала химиче скую систему, т . е . когда значения ф были малы. Однако, как и следовало ожидать, если модель была плохой, проводилось зна
чительно больше итераций, а оценки параметров изменялись |
явно |
случайным образом. Оценки точности вычисленных оценок, |
как |
и оценки точности для Ct, не приводились. |
|
На фиг. П.9.2.1 сравниваются кривые предсказанных концен |
траций, полученные на вычислительной машине, с эксперимен
тальными данными. |
Заметим, |
что основная часть |
наблюдений |
концентраций проводилась |
для установившихся |
состояний |
(t > 1000 мин); такой |
план эксперимента существенно затрудня |
ет оценивание отклика с максимумом, как С2 , или отклика, быстро
|
|
|
|
|
|
|
спадающего |
до нуля, как Су. Тот факт, |
что начальные |
условия |
были известны, существенно помог при оценивании. |
|
|
9.2.3. |
Непрерывные |
наблюдения |
|
|
В примере 9.2.1 было описано оценивание методом наимень |
ших квадратов для случая, когда обработка данных |
проводилась |
с помощью |
цифровой |
вычислительной |
машины, |
не |
связанной |
непосредственно с измерительным устройством. Однако в различ ных отраслях техники, биомедицине, геофизике и океанографии требуется обрабатывать большой объем данных непосредственно в натуральном масштабе времени с тем, чтобы оценивать и регули ровать параметры модели непрерывно, обеспечивая непрерывное уменьшение ф в выражении (9.2.1).
Чтобы непрерывно минимизировать ф, использовалось несколько различных методов, однако здесь будет уделено вни мание лишь наиболее общему из них, а именно методу наиско рейшего спуска. При наискорейшем спуске скорость изменения параметра модели ßj во времени задается соотношением
|
3?---*•&• |
(9 -2 -4) |
Это соотношение выводится следующим образом. |
ф (ß l 5 ß 2 , . . . |
Рассмотрим функцию |
нескольких |
параметров |
. . ., ß n ) . При движении в направлении минимума |
ф необходимо |
пройти в пространстве |
параметров |
некоторое |
расстояние ds, |
