- •1. Основы моделирования динамики систем в программной среде AnyLogic
- •1.1. Общие сведения о возможностях среды AnyLogic
- •1.2. Режим разработки модели
- •1.3. Некоторые специальные приемы разработки моделей в AnyLogic
- •1.4. Режим выполнения модели
- •2. Простые математические модели сложных динамических систем
- •2.1. Модель динамики уровня зеркала водохранилища с гидростанцией
- •2.2. Энергетическая модель сердца
- •2.3. Модель засоления (растворимого загрязнения) водоема с заливом
- •3. Классические эффекты в динамике линейных осцилляторов
- •3.1. Динамика линейного осциллятора
- •3.2. Динамический демпфер
- •3.3. Связанные осцилляторы
- •4. Эффекты в нелинейных динамических системах
- •4.2. Часы Галилея – Гюйгенса
- •4.3. Ламповый генератор электрических колебаний
- •5. Имитационное (агентное) моделирование динамических систем
- •5.2. Модель химической кинетики – брюсселятор
- •Список литературы
В процессе счета (интегрирования) модели по плоскости (u, Q) движется точка, показывающая изменение энергетическо-
го запаса сердца. Под бегунками, изменяющими параметры, находится анимированная картинка сердца, которое бьется при Q > 0 и перестает биться при Q = 0 .
2.2.3. Задания для эксперимента
1. Изменяя параметры функций f (u,Q) и g (u,Q), проследи-
те за изменением графиков этих функций. Сделайте вывод о влиянии этих параметров на характер изменения интенсивностей расхода и пополнения энергетического запаса при изменении управления и самого энергетического запаса.
2.Изменяя значения управления u и начального энергетического запаса Q, проследите за динамикой модели. Как не дать сердцу погибнуть?
3.Изменяя величину a, проследите за изменением «области жизни». Существует ли такое значение a, при котором нормальная работа сердца практически невозможна?
2.3. Модель засоления (растворимого загрязнения) водоема с заливом
2.3.1. Модель
Пусть имеется водоем с заливом, причем вода может перетекать только из водоема в залив, но не в обратную сторону. Примером такой водной системы является Каспийское море с заливом Кара-Богаз-Гол, расположенным примерно на 4,5 м ниже уровня Каспия (рис. 2.6). Пусть V и V1 – объемы воды в море и заливе, Q – интенсивность (объемная скорость) втекания воды в море (реки, дожди, сточные воды), ν – концентрация растворимых примесей (соли и загрязняющих веществ) во втекающем в море потоке воды, I и I1 – интенсивности (объемные скорости) испарения воды из моря и залива, q — интенсивность (объемная
65
скорость) перетекания воды из моря в залив. Кроме того, пусть M и M1 — общие массы растворимого загрязнения в море и заливе, а μ — концентрация загрязнения в море. Считая, что загрязне-
ние в водоеме не достигает насыщения, имеем μ = M /V. В зали-
ве же насыщение имеет место уже давно, и в нем образовались гигантские донные отложения солей натрия, магния и др.
Рис. 2.6
Тогда математическая модель динамики уровней воды и концентраций растворимых примесей в море и заливе имеет вид:
dV |
|
|
|
|
|
|
|
|
= Q − I − q, |
|
|||
dt |
|
|||||
|
|
|
|
|
||
dV1 |
|
|
|
|
||
|
= q − I |
, |
|
|
||
dt |
|
|
||||
|
1 |
|
|
(2.3) |
||
dM |
|
|
|
|
||
|
=νQ − μq, |
|
||||
dt |
|
|
||||
|
|
|
|
|
||
dM1 |
= μq. |
|
|
|
||
|
|
|
||||
|
|
|
||||
dt |
|
|
|
|
|
|
Интенсивность втекания воды в море Q, вообще говоря, может изменяться и зависит от многих причин (изменение температуры, влажности воздуха и др.), но среднее значение Q достаточно постоянно во времени, поэтому будем считать, что Q ≡ const .
Интенсивность испарения воды из моря I является функцией от V, интенсивность испарения воды из залива I1 — функцией от V1, а интенсивность перетекания воды из моря в залив q — функцией
66
от V и V1. Сделаем достаточно очевидное допущение, что функции I (V ) , I1 (V1) и q (V ,V1) обладают следующими свойствами:
I (0) = I1 (0) = 0 , |
dI |
> 0, |
dI1 |
> 0 , |
∂q |
> 0 , |
∂q |
≤ 0 . |
(2.4) |
dV |
|
∂V |
|
||||||
|
|
dV1 |
|
∂V1 |
|||||
По физическому смыслу все фазовые переменные системы (V, V1, M, M1) неотрицательны, поэтому фазовым пространством сис-
темы (2.3) является R 4+ (неотрицательный ортант четырехмерно-
го пространства).
В системе (2.3) первые два уравнения (для V и V1) образуют замкнутую подсистему с фазовым пространством R2+ (первый квадрант плоскости):
|
dV |
|
= Q − I − q, |
|
dV1 |
= q − I1. |
|
(2.5) |
|
dt |
|
|
dt |
|
|||
|
|
|
|
|
|
|
||
Состояния |
равновесия |
(V =V ,V =V ) |
подсистемы (2.5) |
|||||
|
|
|
|
|
|
1 |
1 |
|
ищутся как решения системы уравнений (в общем случае трансцендентных):
I (V ) + q (V ,V1) = Q,
I1 (V1 ) − q (V ,V1 ) = 0,
которая эквивалентна следующей системе (вместо второго урав-
нения напишем сумму обоих уравнений): |
|
|
I(V ) + q(V ,V1 ) = Q, |
(2.6) |
|
I(V ) + I1 (V1 ) = Q. |
||
|
Первое уравнение системы (2.6) определяет некоторую неявную функцию V1 = f1 (V ) , а второе уравнение – неявную функцию
V1 = f2 (V ) .
Покажем, что состояние равновесия подсистемы (2.5) существует и единственно. Дифференцируя второе уравнение (2.6) по V, имеем
dI + dI1 dV1 = 0 , dV dV1 dV
67
и поэтому, в силу (2.4), dV1 / dV < 0 для функции V1 = f1 (V ) , т.е.
эта функция убывающая. Аналогично, из первого уравнения (2.6) имеем
dI + ∂q + ∂q dV1 = 0 , dV ∂V ∂V1 dV
и тогда dV1 / dV > 0 для функции V1 = f2 (V ) , т.е. эта функция возрастающая. Следовательно, графики функций V1 = f1 (V ) и V1 = f2 (V ) могут пересекаться не более чем в одной точке. Будем считать, что f1 (0) ≥ f2 (0) , где V1 = f1 (0) — корень уравнения
I1 (V1 ) = Q , а V1 = f2 (0) — корень уравнения q (0,V1) = Q . Тогда у графиков функций V1 = f1 (V ) и V1 = f2 (V ) существует един-
ственная точка пересечения (V =V ,V1 =V1 ) , а следовательно,
состояние равновесия подсистемы (2.5) существует и единственно.
Теперь покажем, что состояние равновесия (V =V ,V1 =V1 )
глобально асимптотически устойчиво, т.е. при любых начальных условиях из R 2+ траектории подсистемы (2.5) стремятся к состоя-
нию равновесия (V =V ,V1 =V1 ) . Для этого покажем, что: 1) со-
стояние равновесия локально (из малой окрестности) асимптотически устойчиво; 2) в системе не имеется предельных циклов (а также замкнутых контуров, составленных из траекторий).
Локальная асимптотическая устойчивость состояния равновесия проверяется по первому приближению. Сделаем замену пе-
ременных: ξ =V −V , η =V −V ; тогда ξ& |
=V& |
и η& =V& |
, а со- |
|
1 |
1 |
|
1 |
|
стояние равновесия в переменных (ξ,η) имеет координаты (0,0). Далее в новых переменных разложим правые части системы (2.5) в ряд Тейлора по переменным ξ и η в окрестности точки (0,0); получим следующую линеаризованную систему:
68
ξ& = |
|
|
∂q |
|
|
d I |
|
|
|
∂q |
|
|
|
||||
− |
|
− |
|
ξ − |
η, |
|
|||||||||||
|
|
|
|
∂V1 |
|
||||||||||||
|
|
∂V |
|
|
dV |
|
|
|
|
(2.7) |
|||||||
|
|
∂q |
|
|
|
∂q |
|
|
|
dI1 |
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
||||||||
η& = |
|
|
|
|
ξ + |
|
|
|
− |
|
|
η, |
|
|
|||
|
∂V |
|
|
∂V |
|
dV |
|
||||||||||
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
1 |
|
|
1 |
|
|
|
|
||
где производные вычисляются в точке (V =V ,V1 =V1 ) .
Устойчивость или неустойчивость состояния равновесия (0,0) системы (2.7) (а следовательно, и системы (2.5)) определяется
корнями характеристического уравнения λ2 + Aλ + B = 0 . Устойчивость имеет место, если все корни этого уравнения лежат слева от мнимой оси; для этого необходимо и достаточно, чтобы было A > 0 и B > 0 . Вычисляя коэффициенты A и B непосредственно,
имеем (в силу условий (2.4)): |
|
|
|
|
|
|
|
|
|
|
|
||||||||
A = |
∂q |
+ |
|
dI |
− |
|
∂q |
+ |
|
dI1 |
> 0 , |
||||||||
|
|
|
dV |
|
∂V |
|
dV |
||||||||||||
|
|
∂V |
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
1 |
1 |
|
|
|||||
B = − |
dI |
|
∂q |
|
+ |
∂q |
|
dI1 |
|
+ |
dI |
|
dI1 |
> 0 . |
|||||
dV ∂V |
∂V dV |
|
dV dV |
||||||||||||||||
|
|
|
|
|
|
||||||||||||||
1 |
|
|
|
|
|
1 |
|
|
|
|
|
1 |
|
||||||
Следовательно, состояние равновесия локально асимптотически устойчиво.
На плоскости из локальной асимптотической устойчивости состояния равновесия не следует его глобальная устойчивость: в общем случае могут существовать также устойчивые периодические движения. Но в системе (2.5) таких движений нет по критерию Бендиксона. Согласно критерию Бендиксона, если для двумерной динамической системы
V& = F (V ,V1),
V&1 = F1 (V ,V1)
выражение ∂F + ∂F1 сохраняет знак и не обращается тождест-
∂V ∂V1
венно в нуль в некоторой области плоскости (V ,V1) , то в этой
области у системы нет периодических движений (в том числе предельных циклов) и замкнутых контуров, составленных из тра-
69
екторий. Для системы (2.5), в силу сделанных предположений, имеем:
∂F |
+ |
∂F1 |
= − |
d I |
− |
∂q |
+ |
∂q |
− |
d I1 |
< 0 |
|
∂V |
∂V |
dV |
∂V |
∂V |
dV |
|||||||
|
|
|
|
|
|
|||||||
|
|
1 |
|
|
|
|
|
1 |
1 |
|
||
во всём первом квадранте плоскости (V ,V1) . Следовательно, система (2.5) не имеет периодических движений и состояние равно-
весия (V =V ,V =V ) |
глобально асимптотически устойчиво. |
|
1 |
1 |
|
Таким образом, при любых начальных уровнях воды в море и заливе с течением времени обязательно установится водное равновесие.
Процесс установления водного равновесия значительно более быстрый, чем установление равновесной концентрации загрязняющих веществ, поэтому рассмотренная подсистема (2.5) – это не только замкнутая, но и «быстрая» подсистема.
После установления равновесия в «быстрой» подсистеме можно перейти к рассмотрению «медленной» подсистемы (подсистема уравнений для M и M1 в системе (2.3) при условии
V ≡V * и V1 ≡V1* ). Согласно этой подсистеме, равновесная кон-
центрация солей в водоеме μ* = M * / V * удовлетворяет уравне-
нию (условие обращения в нуль правой части уравнения для dM / dt ):
μ* q (V *,V1*) −νQ = 0 ,
откуда находим равновесную концентрацию растворимых солей и примесей в водоеме:
μ* = |
νQ |
|
|
|
. |
(2.8) |
|
q (V *,V * ) |
|||
1 |
|
|
|
Значение μ = μ* является устойчивым состоянием равновесия
для μ , |
так как в силу уравнения для dM / dt в (2.3) величина |
dM / dt |
(а следовательно, и dμ / dt , так как мы рассматриваем |
«медленную» подсистему при условии μ = M / V * ) будет положительной при μ < μ* и отрицательной при μ > μ* .
70
Так как равновесные объемы V * и V1* ищутся как решения системы (2.6), то из (2.8) получаем:
μ* = |
νQ |
|
|
q |
* |
+ |
I |
* |
|
|
I |
* |
|
|
|
I |
* |
|
|
||||
|
=ν |
|
|
=ν 1 |
+ |
|
|
=ν 1 |
+ |
|
|
, |
(2.9) |
||||||||||
|
* |
* |
|
|
|
|
* |
|
|
|
* |
|
|
||||||||||
|
q (V |
) |
|
|
|
q |
|
|
|
|
q |
|
|
|
|
* |
|
|
|||||
|
|
,V1 |
|
|
|
|
|
|
|
|
|
|
I1 |
|
|
||||||||
где звездочки обозначают, что соответствующие величины взяты при равновесных значениях V ≡V * и V1 ≡V1* .
Отношение I * / I1* , входящее в (2.9), грубо приближенно равно отношению площадей водных поверхностей водоема и залива, т.е. S / S1 . Для Каспийского моря и Кара-Богаз-Гола это
отношение примерно равно 37, а следовательно, равновесная концентрация в заливе примерно только в 40 раз больше средней концентрации солей в потоках пресной воды, поступающей в море из рек и с дождями.
Таким образом, залив, подобный Кара-Богаз-Голу, играет роль отстойника, отсасывающего растворимые загрязнения из водоема. Эту же роль в заселяемых стоками воды водоемах могут сыграть искусственные заливы. При этом необходимо обеспечить наличие неисчезающего перепада уровня воды от водоема к заливу, так как обратный поток из залива в водоем сведет на нет отсасывающее действие залива.
2.3.2. Реализация в AnyLogic
Работа реализована в файле Part1\Caspian Sea.alp (рис. 2.7). В окне анимации имеются осциллограммы изменения массы соли и уровня воды в водоеме и заливе. На картинке, схематично изображающей Каспийское море и залив Кара-Богаз- Гол, можно наблюдать подъем и опускание уровней воды. Параметры системы (коэффициенты, входящие в состав функций интенсивностей испарения и притока воды, а также концентрацию соли в речных водах) можно изменять с помощью бегунков.
71
Рис. 2.7
2.3.3. Задания для эксперимента
Изменяя параметры функций интенсивности испарения и притока воды, проследите за изменением динамики системы. Что будет, если искусственно «отрезать» залив от моря с помощью дамбы3 (обращение в нуль коэффициента интенсивности перетекания воды из моря в залив)?
3 Именно так было сделано в 1980-е гг. в отношении Каспия и Кара- Богаз-Гола из опасения, что Каспий обмелеет.
72
