Моделируйте дифференциальные алгебраические уравнения

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

Робертсон [1] создал систему автокаталитических химических реакций для тестирования и сравнения численных решателей для жестких систем. Реакции, константы скорости (k) и скорости реакции (V) для системы даются следующим образом:

$\begin{matrix} A & \overset{k_{1}}{\to} & B & k_{1} = 0.04 & V_{1} = k_{1} [A] \\ B + B & \overset{k_{2}}{\to} & C + B & k_{2} = 3 \cdot 10^{7} & V_{2} = k_{2} [B] [B] \\ B + C & \overset{k_{3}}{\to} & A + C & k_{3} = 1 \cdot 10^{4} & V_{3} = k_{3} [B] [C] \end{matrix}$

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

Simulink Модель из уравнений ОДУ

Система обыкновенных дифференциальных уравнений (ОДУ) имеет следующие характеристики:

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

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

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C^{'} = 3 \cdot 10^{7} B^{2} \end{array}$

Начальные условия: $A = 1$ , $B = 0$ , и $C = 0$ .

Создайте модель

Создайте модель или откройте модель ex_hb1ode.

Добавьте три блока Интегратора к своей модели. Пометьте входы A', B', и C', и выходы A, B, и C соответственно.
Добавьте блоки Sum, Product и Gain, чтобы решить каждую дифференциальную переменную. Для примера, чтобы смоделировать сигнал C',
1. Добавьте блок Math Function и соедините вход с сигналом B. Установите параметр Function равным square.
2. Соедините выход из блока Math Function с блоком Gain. Установите параметр Gain равным 3e7.
3. Продолжите добавлять оставшиеся условия дифференциального уравнения к вашей модели.
Моделируйте начальное условие A путем установки параметра Initial condition для A Блок интегратора для 1.
Добавьте блоки Out, чтобы сохранить сигналы A, B, и C в переменный MATLAB yout.

Симулируйте модель

Создайте скрипт, который использует sim команда для симуляции модели. Этот скрипт сохраняет результаты симуляции в переменном MATLAB yout. Поскольку симуляция имеет длительный временной интервал и B первоначально изменяется очень быстро, графическое изображение графиков значений в логарифмической шкале помогает визуально сравнить результаты. Кроме того, поскольку значение B является маленьким относительно значений A и C, умножить B около $1 \cdot 10^{4}$ перед построением графика значений.

Введите следующие операторы в MATLAB^® скрипт. Если вы создали свою собственную модель, замените ex_hblode с именем вашей модели.
```
sim('ex_hb1ode')
yout(:,2) = 1e4*yout(:,2);
figure;
semilogx(tout,yout);
xlabel('Time');
ylabel('Concentration');
title('Robertson Reactions Modeled with ODEs')
```
Из Simulink^® Редактор, на вкладке Modeling, нажмите Model Settings.
- На панели Solver установите значение Stop time 4e5 и Solver для ode15s (stiff/NDF).
- На панели «Импорт данных» установите флажки Time и Output.
Запустите скрипт. Заметьте, что все A преобразуется в C.

Simulink Модель из уравнений ДАУ

Система дифференциальных алгебраических уравнений (ДАУ) имеет следующие характеристики:

Он содержит как обыкновенные дифференциальные уравнения, так и алгебраические уравнения. Алгебраические уравнения не имеют никаких производных.
Только некоторые уравнения являются дифференциальными уравнениями, определяющими производные некоторых зависимых переменных. Другие зависимые переменные заданы алгебраическими уравнениями.
Количество уравнений равно количеству зависимых переменных в системе.

Некоторые системы содержат ограничения из-за законов сохранения, таких как сохранение массы и энергии. Если вы устанавливаете начальные концентрации равными $A = 1$ , $B = 0$ , и $C = 0$ , общая концентрация трех видов всегда равна $1$ с тех пор $A + B + C = 1$ . Вы можете заменить дифференциальное уравнение для $C^{'}$ с помощью следующего алгебраического уравнения для создания множества дифференциальных алгебраических уравнений (ДАУ).

$C = 1 - A - B$

Дифференциальные переменные A и B однозначно определите алгебраическую переменную C.

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C = 1 - A - B \end{array}$

Начальные условия: $A = 1$ и $B = 0$ .

Создайте модель

Внесите эти изменения в свою модель или в модель ex_hb1ode, или откройте модель ex_hb1dae.

Удалите блок Интегратор для вычисления C.
Добавьте блок Sum и установите параметр List of signs равным + - -.
Соедините сигналы A и B к минусовым входам блока Sum.
Моделируйте начальную концентрацию A с блоком Constant, соединенным с входом плюс блока Sum. Установите параметр Constant value равным 1.
Соедините выход блока Sum с разветвленной линией, соединенной с блоками Product и Out.

Симулируйте модель

Создайте скрипт, который использует sim команда для симуляции модели.

Введите следующие операторы в скрипт MATLAB. Если вы создали свою собственную модель, замените ex_hbldae с именем вашей модели.
```
sim('ex_hb1dae')
yout(:,2) = 1e4*yout(:,2);
figure;
semilogx(tout,yout);
xlabel('Time');
ylabel('Concentration');
title('Robertson Reactions Modeled with DAEs')
```
В редакторе Simulink Editor на вкладке Modeling нажмите Model Settings.
- На панели Solver установите значение Stop time 4e5 и Solver для ode15s (stiff/NDF).
- На панели «Импорт данных» установите флажки Time и Output.
Запустите скрипт. Результаты симуляции, когда вы используете алгебраическое уравнение, те же, что и для симуляции модели, используя только дифференциальные уравнения.

Simulink Модель из уравнений ДАУ, использующих Алгебраический Ограничительный Блок

Вы можете заменить дифференциальное уравнение для $C^{'}$ с алгебраическим уравнением, смоделированным с помощью блока Алгебраических Ограничений и блока Суммы. Блок Алгебраического Ограничения ограничивает свой входной сигнал F (z) на нуль и выводит алгебраическое состояние z. Другими словами, выход блока является значением, необходимым для получения нуля на входе. Используйте следующее алгебраическое уравнение для входа в блок.

$0 = A + B + C - 1$

Дифференциальные переменные A и B однозначно определите алгебраическую переменную C.

$\begin{array}{l} A^{'} = - 0.04 A + 1 \cdot 10^{4} B C \\ B^{'} = 0.04 A - 1 \cdot 10^{4} B C - 3 \cdot 10^{7} B^{2} \\ C = 1 - A - B \end{array}$

Начальные условия: $A = 1$ , $B = 0$ , и $C = 1 \cdot 10^{- 3}$ .

Создайте модель

Внесите эти изменения в свою модель или в модель ex_hb1ode, или откройте модель ex_hb1dae_acb.

Удалите блок Интегратор для вычисления C.
Добавьте блок Алгебраических ограничений. Установите параметр Initial guess равным 1e-3.
Добавьте блок Sum. Установите параметр List of signs равным - + + +.
Соедините сигналы A и B в плюс входы блока Sum.
Моделируйте начальную концентрацию A с блоком Constant, соединенным с минусовым входом блока Sum. Установите параметр Constant value равным 1.
Соедините выход блока Алгебраических ограничений с разветвленной линией, соединенной с входами блоков Product и Out.
Создайте линию ветви от выхода блока Алгебраических ограничений до последнего плюс входа блока Sum.

Симулируйте модель

Создайте скрипт, который использует sim команда для симуляции модели.

Введите следующие операторы в скрипт MATLAB. Если вы создали свою собственную модель, замените ex_hbl_acb с именем вашей модели.
```
sim('ex_hb1dae_acb')
yout(:,2) = 1e4*yout(:,2);
figure;
semilogx(tout,yout);
xlabel('Time');
ylabel('Concentration');
title('Robertson Reactions Modeled with DAEs and Algebraic Constraint Block')
```
В редакторе Simulink Editor на вкладке Modeling нажмите Model Settings.
- На панели Solver установите значение Stop time 4e5 и Solver для ode15s (stiff/NDF).
- На панели «Импорт данных» установите флажки Time и Output.
Запустите скрипт. Результаты симуляции, когда вы используете блок Алгебраических ограничений, те же, что и для симуляции модели, используя только дифференциальные уравнения.

Использование блока Алгебраического Ограничения создает алгебраический цикл в модели, Если задать параметр Algebraic Loop равным warning (на вкладке Modeling нажмите Model Settings, затем выберите Diagnostics), следующее сообщение отобразится в Средство Просмотра во время симуляции.

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

Ссылки

[1] Робертсон, H. H. «решение набора уравнений скорости реакции». Численный анализ: введение (J. Walsh ed.). Лондон, Англия: Академическая пресса, 1966, стр. 178-182.

Документация