CC BY
139
УДК 519.6
Я.А. Холодов, А.С. Холодов, А.В. Гасников, И.И. Морозов, В.Н. Тарасов
Московский физико-технический институт (государственный университет)
Моделирование транспортных потоков — актуальные проблемы и перспективы их решения
Работа посвящена математическому моделированию транспортных потоков, описанию проблем, которые при этом возникают, и возможностей их решения. Предложена оригинальная макроскопическая модель для описания автомобильного движения на сложном графе транспортной сети. Макроскопическая модель основана на гидродинамическом подходе, где транспортный поток описывается уравнениями течения сжимаемой многокомпонентной жидкости с мотивацией. В качестве компоненты рассматриваются автомобили, объединённые общими характеристиками. Разработан также оригинальный алгоритм построение систем уравнений в узлах графа транспортной сети — перекрестках. Результаты численных расчётов при сравнении с экспериментальными данными показывают работоспособность предлагаемой модели.
Ключевые слова: автомобильное движение, гидродинамический подход, граф транспортной сети, многокомпонентные транспортные потоки.
I. Введение
Транспортные проблемы современных мегаполисов хорошо известны и достаточно актуальны. Решение этих проблем не всегда очевидно и уж точно всегда требует времени, при этом никто не гарантирует, что в процессе решения текущих проблем не возникнут новые. В этой связи возникает вопрос: как можно облегчить жизнь непосредственным участникам дорожного движения, а именно водителям. Ответ на этот вопрос не однозначен, поскольку требует комплексного подхода и глубокого понимания существующих на сегодняшний день проблем. Реально складывающаяся дорожно-транспортная ситуация требует комплексных решений не только в сфере развития и оптимизации транспортной инфраструктуры, но и в повседневном управлении транспортной системой для обеспечения надежности её функционирования и безопасности её участников как в штатных, так и в кризисных ситуациях. При этом важно учитывать не только текущую дорожную обстановку, но и иметь детальный ситуативный прогноз на всю транспортную систему с учетом временных и количественных затрат реализации различных сценариев развития. Для этого необходима корректная вычислительная математическая модель, которая будет осуществлять данный прогноз по всей транспортной сети мегаполиса.
Прежде всего следует определить модель и её параметры, с которыми предстоит работать. В разрабатываемой модели следует учесть все возможные факторы, такие как периодически обновляемые данные о текущем состоянии дорожной ситуации, правила движения по полосам на перекрестках, расписания и алгоритмы работы светофоров, локальные перекрытия на участках дороги и т.д. Модель также должна включать временные
зависимости коэффициентов распределения потоков на перекрестках — матрицы перемешивания вместе с матрицей корреспонденций для заданной транспортной сети.
Транспортный поток может моделироваться в зависимости от плотности и скорости дорожного движения разными способами. Начиная от простейших макроскопических моделей, используемых в основном для однополосного беззатор-ного движения автомобильного транспорта, и заканчивая моделированием движения в многополосной пробке на основе гидродинамического приближения, когда транспортный поток уподобляется сжимаемой многокомпонентной жидкости с мотивацией, учитываемой через уравнение состояния транспортного потока. Ключевым моментом также является алгоритм построение систем уравнений в узлах графа транспортной сети — перекрестках. Данные системы уравнений обеспечивают связь моделируемых величин на всем графе транспортной сети путем задания корректных граничных условий для всех его ветвей, входящих и выходящих в каждый из его узлов.
Отдельно встает задача о выборе оптимального маршрута для каждого из водителей, которую можно рассматривать как задачу принятия решения. А именно, путем исследования множества достижимости куда и с какими характеристиками можно доехать, например, из заданной точки за данное время. При этом важно учитывать не только текущую дорожную обстановку, но и иметь детальный прогноз на всю траекторию движения с учетом временных затрат по маршруту. Мысленно себе это можно представлять как исследование растекания жидкости из вершины по транспортному графу.
В основу предлагаемой технологии положен имеющийся опыт использования сетевых вычис-
лительных моделей, основанный на решении соответствующих краевых задач для уравнений в частных производных на графах [1, 2]. Данный опыт показывает, что при современном уровне развития вычислительной математики и вычислительной техники, особенно с появлением высокопроизводительных вычислительных платформ и решений, такой подход позволяет эффективно решать достаточно сложные «глобальные» задачи данного класса [3].
II. Сетевая модель интенсивного дорожного движения в мегаполисе
Для макромоделирования автомобильного движения используются различные математические модели, в том числе основанные на дифференциальных уравнениях в частных производных и, в частности, гидродинамические модели, аналогичные уравнениям течения сжимаемой многокомпонентной жидкости с мотивацией (см. [4-11] и др.). Ниже также рассматривается вычислительная математическая модель интенсивного уличного движения в мегаполисе, основанная на решении соответствующих краевых задач для уравнений в частных производных гиперболического типа.
Будем представлять уличную сеть мегаполиса в виде перекрестков (узлов графа с номерами I — 1, Ь). связанных между собой доро-
гами с длиной Хк (ребрами графа с номерами к — 1, ..., К, 0 ^ Хк ^ Хк) (см., например, рис. 1). Введем следующие обозначения для используемых в модели параметров: рк (Хк,і) — плотность автомобильного потока (количество автомобилей на единицу длины дороги с номером к);
Vк (Хк ,£) — средняя по сечению дороги линейная скорость потока автомобилей на к-й ветви графа в точке Хк в момент времени £; дк (Хк ,ї) — РкVк — поток автомобилей.
l=L+L к=10 l=L+L+3... к=8
^; ,*fX7 х=0.
L+LV к=к l=L к=9 1=3... к=7 к=6 l=L+3
l=L+1 к=11...(т=2) 1=1 1=2
к=1(т=1) к=3 (т=3)... к=2 (т=М) l=L+£+1 к=5 к=4 l=L+2 l=L+L+2
Рис. 1. Пример направленного графа для перекрестка
Н.1. Система уравнений автомобильного движения на ребре графа транспортной сети
Гиперболическая система уравнений, описывающая автомобильное движение, представляет со-
бой дифференциальные законы сохранения (изменения) «массы» и «импульса» на автодорогах (по аналогии с гидродинамикой) и в дивергентной форме записывается в виде (см., например, [6]):
Г dpk/dt + д(ркVk)/дхк = fok, m
\ д(ркVk)/dt + д(ркvk + p(pk))/дхк = f\k. ( 1
Здесь
p = p(pk1 () есть замыкающее систему (1) уравнение состояния (зависимость давления от плотности, рис. 2), fok — возможные источники или стоки «массы» (въезжающие на дорогу автомобили из не учитываемых явно элементов уличной сети, останавливающиеся или начинающие движение автомобили и т.п.), fik — возможные внешние импульсы, действующие на систему (1). В таких моделях основной проблемой является построение адекватного действительности уравнения состояния (2), конкретный вид которого, как и для всякой феноменологической модели, должен быть определен из экспериментальных измерений (возможно с использованием параметрических решений системы (1), (2)). Представленные на рис. 3, 4, 5, 6 для движения на одной (левой) полосе высокоскоростной автострады экспериментальные зависимости потока от плотности из [12-15] показывают (отмечены точками и соединяющими их отрезками прямых), что предельная величина потока (определяемая пропускной способностью дороги, техническими характеристиками автомобилей и правилами дорожного движения) много больше реально реализуемой максимальной его величи-
(free)
ны qfs = Qmax и реально реализуемое движение в фазовой плоскости {ps,qs} (номер ветви графа к далее опущен) имеет своей границей параболу (кривая № 1 на рис. 3-6):
qs(ps1 = ps(qfs/pfs + Cfs(l - Ps/Pfs)),
Ps < Pfs (3)
при плотности одной полосы дороги, меньшей критической ps < pfs (область «свободного потока»). Здесь Cfs — аналог скорости звука при критической плотности, s — номер полосы движения (1 ^ s ^ S), S — число полос заданной автодороги.
0 0.2 0.4 0.6 0.8 р/р.
Рис. 2. Примеры уравнения состояния р(р)
q
0 20 40 60 Ш р ЦНУЩ
r ' km '
Рис. З. Аппроксимация экспериментальных данных из [1З]
І 20 40 60 p^iggssj
Рис. 4. Аппроксимация экспериментальных данных из [1З]
q vehicles li
1 20 40 SO 80 рГ^
Рис. 5. Аппроксимация экспериментальных данных из [14]
Однако при дальнейшем увеличении плотности до значений ps > pfs такой однозначной зависимости, как в (3), из представленных экспериментальных данных не наблюдается вследствие взаимодействия волн разгона и торможения при интенсивном (с закритической плотностью ps > pfs) автомобильном движении. Скорость свободного движения потока автомобилей vs0(0) = (dqs(ps)/dps)0 определяется правилами дорожного движения, техническими характеристиками автомобилей и стилем вождения участников дорожного движения и обеспечивает связь этого параметра со значениями qfs,pfs,Cfs:
Vso = Vs(0) = qfs/pfs + Cfs. (4)
Предполагая существование аналогичной (3) однозначной зависимости потока от плотности также и при ps > pfs, будем полагать ее линейной в
области реализуемых значений при ps > pfs (прямая № 2 на рис. 3-6):
qs(ps) = C*s(P*s Ps^l qs(pfs) = C*s(P*s Pfs) =
-- dsqfs, Ps > Pfs (5)
с обращающимся в ноль при максимальной плотности ps = p*s потоком и с возможным его разрывом в случае ds < 1 при критической плотности ps = pfs, коэффициент разрыва ds определяется из измеренных данных:
ds — C*s (P*s Pfs)/qfs ^ 1■ (6)
Здесь p*s — предельно допустимая плотность по-
тока автомобилей (пропускная способность одной полосы дороги), соответственно C*s — аналог скорости звука при этой плотности.
q (vechicles/li) DI-left lane
0 10 20 30 40 50 p (vechicles/km)
Рис. б. Аппроксимация экспериментальных данных из [14]
Зависимости qs = qs(ps) из (З), (5) или аналогичных аппроксимаций измеренных данных вместе с первым уравнением из (1) (при fok = 0 дифференциальным законом сохранения числа автомобилей), так же как и система уравнений (1)-(2), могут непосредственно использоваться для моделирования дорожного движения. В этом случае движение автомобилей описывается одним нелинейным уравнением переноса для плотности (см., например, [4, 5]):
dps/dt + dqs(ps)/ дх = dps /dt+ Щ^-др3/дх =
dps
= 0. (7)
Скорость Vs всегда может быть вычислена с использованием рассчитанных из (7) значений ps и известной фундаментальной диаграммы qs(ps):
Vs(x,t) = qs (ps (x^^s^^). (S)
В данной работе ставится тем не менее другая задача — получить из известного вида фундаментальной диаграммы аналог уравнения состояния
(зависимость давления от плотности) для заданной полосы автодорогира = ра(ра) и использовать для решения системы (1). Использование одного уравнения не является достаточным для корректного описания всех фаз (состояний) транспортного потока, наблюдаемых на рис. 3-6 (см. [11]). Сделать это можно, если воспользоваться хорошо известным дифференциальными преобразованиями законов сохранения (см., например, [16]). Умножим (7) на дда/дра:
^-др а/д1 + ^-дда/дх = дда/сП. + ^-дда/дх =
дра дрэ дрэ
0
(9)
и воспользуемся следующими дифференциальными соотношениями:
дда/дЬ = д(реуе)/дг = Уа дра/дЬ + радУа/дЬ,
дд,э/дх = д(р8у8 )/дх = Уа дра/дх + р8ду8/дх, подставим их в (9), откуда получим
ду3/д1 -|—Щ^^-дг13/дх = 0.
дра
(10)
Итак, вместо системы уравнений (1) получен аналог этих законов сохранения в виде системы уравнений (7), (10) в недивергентной форме с нулевой правой частью:
Г дрз/дг+ЩеАдрз/дх = о,
\ ду3/д1+Щ^ду3/дх = 0. 1 ]
Из сопоставления уравнений (10) и второго уравнения в (1), записанного в недивергентной форме, можно получить вид уравнения состояния (2). Для этого в (10) добавим и вычтем член уадуа/дх и воспользуемся следствием второго дифференциального соотношения: дУа/дх = (дда/дра — да/ра)(дра/дх)/ра. Тогда (10) приводится к виду:
дуа /дЬ + Уа дуа/дх +
(да/ра - дда/дра)2 Ра
дра/дх = 0.
Сравним это уравнение с недивергентной формой второго уравнения системы (1) для одной полосы с нулевой правой частью: дуа/д1: + г>адг>а/дх + -^-щ-дра/дх = 0, получим соотношение, связывающее производную, давление дра/дра, аналог скорости звука са, поток да и плотность ра:
дра/дра = с23 = (да/ра — дда/дра)2. (12)
Интегрируя (12), можно найти функциюра(ра) по заданной функции да(ра), если последняя описывает достаточно широкий класс дорожного движения.
Для да(ра) из (3), (5) имеем
Ра (ра)= с^ар3а/(3р^а) при ра < Г^а, (13)
Ра(ра) = Р!а + р* ас* а(1/р !а — 1/Pа),
Р!а = с)ар!а/3 при ра >р/а. (14)
Здесь и далее мы полагаем вектор правой части ] = 0 из следующих соображений. Зависимости потока от плотности, показанные на рис. 3-6, получены от датчиков, которые находились в левой полосе автобана на существенном расстоянии от каких-либо дорожных развязок (см. [13-14]), соответственно там не было никаких возможных источников или стоков «массы» автомобилей. Все остальные параметры, которые исследователи обычно так любят включать в правую часть системы уравнений (1) (см. [4-10] и др.), учитываются автоматически через полученный вид фундаментальной диаграммы для этих измерений и в правую часть не входят.
Вблизи предельных значений плотности ра ~ р*а движение является слабосжимаемым, поэтому для обеспечения бесстолкновительного движения в области сильных заторов бралось р*а ^ ра ^ ре = р*а(1 — е), £ ~ 0,1 — 0,01, соответственно уравнение состояния (14) модифицировалось следующим образом:
ра(ра) = (Реа + (ра — Г->еа )(В + Ср + реа))),
С = ((др/др)*а — (др/др)еа)/(р*а — ре),
В = (др/др)*а — 2Ср*.
(15)
Выбором значений ре и (др / др)=ас2Еа можно также обеспечить приемлемый шаг интегрирования по времени в области автодороги, где плотность приближается к своему предельному значению р*а, а значение (др/др)*а начинает неограниченно расти.
Чтобы перейти в системе уравнений (1)-(2) от многополосного к усредненному представлению дороги, для ветви графа с номером к = 1, ..., К в (1)-(2) аналогично многокомпонентной односкоростной модели баротропного газа будем полагать
1
вк
г’к = Yk^Vsk,
а= 1
(16)
вк вк
дк = &кдк ''У ^ дак , рк = &крк ^ ^ рак .
а=1 а=1
Переход (16) удобен тем, что его можно осуществлять как в прямую, так и в обратную сторону в любом участке дороги при возникновении такой необходимости. Индекс в можно связывать не только с номером полосы дороги, но и с маршрутом водителя, типом автомобиля, стилем вождения или опытностью водителей и т.д.
При этом для различных участников движения, имеющих существенно отличающиеся друг от друга константы в уравнениях состояния (13),
(14), в дополнение к системе (1), (2), как это принято, можно ввести в рассмотрение их концентрации:
вк
Сак = рак /'^ рак, а=1
вк
Е
а=1
Сак = 1
(17)
в такой равновесной или неравновесной, односкорОСТНОЙ (V/. = г’зк) ИЛИ многоскорост-
ной смеси и в простейшем случае для определения этих концентраций использовать соответствующие уравнения конвективно-диффузионного переноса (в том числе между разными полосами дороги):
д2Сак
+ /(С1к, ..., Свк), (18)
дСзк , дСзк
сН. дх А‘ дх1
где р — «коэффициент диффузии», а /(С1к, ..., Свк) определяет, например, процессы перестроений и т.п.
Имея замкнутую систему уравнений (1)--(2) (хорошо известную систему уравнений Эйлера для сжимаемого баротропного газа), можно найти уравнения характеристик и соответствующие условия совместности этой системы уравнений гиперболического типа, соотношения на разрывах, возникающих при переходе от докритического («сверхзвукового») движения к закритическому («дозвуковому») в волнах торможения и т.д. Интегрируя (1), (2) по замкнутому контуру в плоскости {х,Ь}, содержащему фронт разрыва движущегося со скоростью Б, можно получить соотношения на разрывах, в частности, уравнение ударной адиабаты, связывающее параметры перед фронтом разрыва (в докритическом потоке с р1к ^ р/к, кривая № 1 на рис. 3, 4, 5, 6, 7, 8) и за фронтом
разрыва (в закритическом потоке с р2к > р/к, прямая № 2 на рис. 3-8):
(Р1к (р1к ) — Р2к (р2к))(1/р2к — 1/р1к) — (У1к — У2к) = = 0, (19)
а также скорость движения фронта разрыва:
Б = (р1кУ1к — р2кУ2к )/(р1к — р2к) =
= (д1к — д2к)/(р1к — р2к). (20)
Эти соотношения и, в частности, вполне измеряемую величину Б можно использовать для определения констант в уравнении состояния (13)--(14), как это делается при построении уравнений состояния сложных физических сред.
Рис. 7. е/в = 0.1[м/с]: предельные зависимости потока от плотности — кривые № 1, 2; ударные адиабаты — сплошные кривые без номеров; характеристики первого семейства — штриховые линии
Рис. 8. = 0.5: предельные зависимости потока от
плотности — кривые № 1, 2; ударные адиабаты — сплошные кривые без номеров; характеристики первого семейства — штриховые линии
Вводя в рассмотрение векторы V = {ркУк},
и(V) = {рк,ркУк}, Г(V) = {ркУк,рк+ р(рк)},
/ = {/ок,/1к}, систему уравнений (1) можно записать в векторной дивергентной форме:
ди^ )/дЬ + дГ (V )/дхк = / (У,хк,Ь), или в недивергентной форме:
ЗУ/дЬ + А(V )дУ/дхк = / = (дU/дV )-1 / с матрицей Якоби:
А = (дU/дV )-1(дF/дV) =
(21)
Ук рк
(др/дрк )/рк Ук
(22)
неособенная, так
Собственные числа матрицы Якоби А : Ах = г’к+Ск, А2 = 1’к - ск, ск = Vхдр/дрк всегда действительные и различные (др/дрк = с\ > 0) —
система (21) строго гиперболического типа, а левые собственные векторы: ш1 = {ск,рк},
ш2 = {ск, — рк} линейно независимы. Т. е. матрица со строками из левых собственных векторов
ш1 \ / ск рк
Ш2 ) \ ск —рк
как ёе! П = —2ркск = 0 при ненулевых плотности рк и «скорости звука» ск. Для упрощения выкладок в дальнейшем будем полагать вектор правой части / = 0.
Эквивалентная (1) характеристическая форма уравнений совместности вдоль характеристик первого (верхний знак, индекс № 1 ) и второго (нижний знак, индекс № 2) семейств
дх = \дЬ = (ук ± ск)дЬ, г = 1, 2
имеет вид
Ск(,Рк)дрк/д£-1 гЬ ркдУк/д£-1 — — — + А-— г — 1 9
<1и ~ дг ^ л> дхк > * —
(23)
Каждое из уравнений (24) является, по сути, обыкновенным дифференциальным уравнением вдоль i-й характеристики (23). Для уравнения состояния (1З)-(14) условия совместности (24) интегрируются аналитически, что позволяет получить вид инвариантов Римана ri, постоянных вдоль соответствующих характеристик Xi, i = І, З:
n(pk,Vk) = CfsPk|(Spfs) ± Vk = const,
при pk < Spfs; (ЗБ)
ri(Pk,Vk) = ~(c*sP*sS/pk) ± Vk = const,
при Pk > Spfs■ (26)
В фазовой плоскости {p,q} кривые (25), (26) имеют вид
qk(pk) = ±Pk(ri - CfsPk/(Spfs)) = const, i = 1.2 при pk < Spfs; (27)
qk(pk) = ±(npk + C*sp*sS) = const, i = 1.2
при pk > Spfs. (28)
Характеристическая форма уравнений (24) или ее модификации с использованием инвариантов Римана (25), (26) (если удается найти точное решение (24)), так же как соотношения на разрывах (19), (20), необходимы при построении эффективных разностных схем для численного решения системы (1)-(2) (см., например, [17]), а в граничных точках они используются для замыкания краевых условий.
Для представленных на рисунках экспериментальных данных из [13] (рис. 3—4) и из [14] (рис. 5-6) (точки и соединяющие их отрезки прямых) константы pfs,qfs,Cfs,ds,p*s,c*s в уравнениях (3), (5), аппроксимирующих границы области реализуемых движений автомобилей, принимались следующими.
На рис. 3:
pfs = 0,0243[АТС/м],
qfs = 0,656[АТС/с], Cfs = 15,8[м/с],
< ds = 1, (29)
p*s = 0,21[АТС/м], Cts = 3,53[м/с], vso = 42,8[м/с].
На рис. 4:
pfs = 0,0243[АТС/м],
qfs = 0,656[АТС/с], Cfs = 27[м/с], ds = 1, p*s = 0,21[АТС/м], c^s = 3,53[м/с],
Vs0 = 54[м/с].
На рис. 5:
Р}3 — 0,0213[АТС/м],
qfs — 0,6[АТС/с], си — 22,4[м/с],
< ds — 0,9, (31) рts — 0,21[АТС/м], cts — 2,9[м/с],
Vs0 — 50,6[м/с].
На рис. 6:
р]:.э — 0,0213[АТС/м],
qfs — 0,6[АТС/с], cfs — 22,4[м/с],
< ds — 0,9, (32) р^s — 0,21[АТС/м], c■ts — 3,17[м/с],
Vs0 — 50,6[м/с].
Соответствующие им граничные зависимости потока от плотности (3), (5) показаны сплошными линиями 1, 2 на рис. 3-8. Для этих же значений констант уравнения состояния (13), (14) представлены на рис. 2. Кривые а, Ь, с, d соответствуют константам из (29), (30), (31) и (32).
На рис. 3-6 сплошными линиями без номеров показаны ударные адиабаты (19) для ряда значений р\к^\к — дік/рік, принадлежащих кривым (3), а штриховыми линиями показаны инварианты Римана первого семейства (27), (28), і — 1. Инварианты Римана второго семейства і — 2 симметричны первому семейству относительно оси абсцисс. Отдельно стоит обратить внимание на рис. 5. Здесь мы видим, что переход из докритиче-ского («сверхзвукового», кривая № 1) движения к закритическому («дозвуковому», прямая № 2) осуществляется через ударную волну сжатия и коэффициент разрыва дц < 1. Таким образом, данный процесс не является адиабатическим. Это приводит к тому, при возвращении назад в ту же точку на кривую № 1 с прямой № 2 через волну разрежения невозможно.
Для оценки влияния констант в уравнении состояния (13), (14) на поведение давления в варианте (30) более существенно, чем в (29), (31) и (32), были изменены коэффициенты cfs,ds. На рис. 7 для случая cfs — 0,1 [м/с] и на рис. 8 для
— 0,5 представлены аналогичные описанным выше зависимости параметров: предельные зависимости потока от плотности — кривые № 1, 2; ударные адиабаты — сплошные кривые без номеров; характеристики первого семейства — штриховые линии. Соответствующие зависимости для давления показаны на рис. 2 (кривые е и ]). Видно, что существенное отличие параметра от 1 приводит к заведомо неправильному поведению ударных адиабат (увеличение скорости при переходе от докритического режима к закритическому в волне торможения и т.п.). Из сравнения кривых а и е на рис. 2 видно также, что аппроксимация в (3) докритической ветви дДрл) даже линейной функцией (cfs « 0) влияет на уравнение состояния сравнительно слабо. Сопоставление кривых а и c на рис. 2, для которых дqs(рfs)/дрs > 0,
с кривыми Ь и д, для которых дда (р/а)/дра = 0, показывает, что равенство нулю дда(р/а)/дра = 0 дает более гладкое сопряжение докритической и закритической ветвей ра(ра) и обеспечивает принципиальную возможность перехода в закритиче-ский режим движения без экстренного торможения (аналог адиабатического безударного сжатия потока в газовой динамике). В этом случае закри-тическая ветвь зависимости потока от плотности (5) (прямая № 2 на рис. 6) является слабой волной торможения (19) (при переходе через которую параметры непрерывны, а разрывными являются нормальные производные) и при этом остается одной из характеристик семейства № 1.
Н.2. Система уравнений автомобильного движения в узлах графа транспортной сети
Число граничных условий на свободных концах ветвей графа дорожной сети (входах хк = 0 и выходах хк = Хк) зависит от знаков собственных чисел матрицы Якоби А (22). Их количество на входах может быть равным двум при положительных А1 ,\2 > 0, что полностью определяет параметры такого узла, одному при А1 > 0,А2 < 0 или нулю при отрицательных А1Д2 < 0. На выходах их число может быть равным нулю при положительных А1Д2 > 0, единице при А1 > 0,А2 < 0 или двум при отрицательных А1А < 0.
В соответствии с этим на входах-выходах из дорожной сети в качестве граничных условий могут быть заданы как функции времени значение плотности потока автомобилей рк(Ь) или величина линейной скорости Ук (Ь) (при А1 > 0,А2 < 0), для расчета таких узлов в качестве второго уравнения привлекается одно из условий совместности
(24) или (25), (26) вдоль идущих внутрь области интегрирования характеристик. Если необходимо задавать два граничных условия, обе эти переменные задаются одновременно. Также возможна ситуация, когда граничные узлы рассчитываются через значения во внутренних точках дороги, тогда задавать граничные условия не требуется. Вместо плотности или линейной скорости может задаваться величина автомобильного потока
дк (Ь) = ркик . (33)
На выходах иногда используют неотражающие граничные условия нулевые производные
(ддк/дхк) = 0, (дик/дхк) = 0, (34)
что, вообще говоря, является корректным лишь при положительных А1Д2 > 0.
Помимо граничных условий, для системы (1), (2) необходимо задать также некоторые начальные условия:
В точках ветвления графа (во внутренних узлах или на перекрестках) постановка граничных условий существенно зависит от организации движения (учитывается ли многополостное движение на дороге или рассматривается осредненное по направлению движение, учитывается ли наличие светофоров или рассматриваются более длительные временные масштабы и т.д.).
В каждом внутреннем узле графа дополнительной искомой величиной, помимо значений векторов V11, V12, ..., Vim в окончаниях входящих в узел l и выходящих из него m = 1, ..., M ветвей, является «давление» узла pi(t). При расчете таких граничных точек в зависимости от знаков собственных чисел Ai,A2 матрицы A для различных связанных с узлом l ветвей m = 1, ..., M возможны различные ситуации. При Ai > 0,А2 < 0, то есть при \vik \ < cik (закритическое движение), в дополнение к условиям совместности (24) (или
(25), (26)) вдоль направленных внутрь области интегрирования характеристик (23) в качестве граничных условий использовались соотношения
Ф1т (t,pl,Vlm ) = plm pl Rlmplmvlm = 0. (36)
Здесь Rlm ^ 0 — коэффициент сопротивления окончания m-й ветви, входящей в l узел (или выходящей из него), plm = p(plm) — «давление» в примыкающем к узлу l конце ветви m, связанное с плотностью plm в этой точке уравнением состояния (2) (или (13), (14)). С помощью коэффициентов Rlm(t,plm,Vlm), в частности, можно воспроизводить правила проезда нерегулируемых перекрестков, работу светофоров и т.п.
Помимо аналогов закона Ома (36) или аналогичной формулы Пуазейля для течений в трубах, для расчета внутренних узлов графа используется равенство нулю алгебраической суммы потоков автомобилей:
м
^l (t,Vl1, ..., VlM ) ^ ^ Plm plmvlm = °- (37)
т=1
В (37) коэффициенты f3lm = 1 для входящих в узел ветвей и f3lm = —1 для выходящих из узла ветвей.
Таким образом, M дифференциальных условий совместности (24) (или ((25), (26)) и M +1 нелинейных алгебраических граничных условий (36), (37) позволяют (после их аппроксимации разностными соотношениями) найти 2M +1 искомых параметров Vn = {pn,vn}, ..., Vim = {pimVim}, pl (t) . Для решения получившейся нелинейной системы уравнений можно воспользоваться каким-либо итерационным методом. Например, аппроксимируя (24) (или ((25), (26)) по какой-либо явной схеме:
<Plm (Vlm ) = ^Im (V^ — V^m ) = cU^ — Pkm)±
Pk (xk ,0) = Pko(xk), Vk(xk, 0) = Vko (xk). (35)
±Pkm(vk+1 — vkm) =0, m =1, ..., M. (38)
В (38) У^т является значением сеточной вектор-функции в точке пересечения характеристики (23), приходящей в рассчитываемую точку временного слоя Ь = 1п+1, со слоем Ь = 1п для к-го ребра. Эта величина определяется интерполяцией по данным вектор-функции УП в узлах сетки на слое £ = Ьп. В целом, вместо (23) предпочтительнее использовать разностную аппроксимацию соответствующего инварианта Римана (25), (26), обеспечивающую более высокую точность расчетов.
Обозначив вектор неизвестных как УI = {Р11, VII, ■ ■ ■, Р1М , У1М , Рь}, получаем систему нелинейных уравнений:
(У) = {^1 (Уц, ■ ■ ■, Ум), фп(Уп), ■ ■ ■, фхм(Ум),
Фи(рь Уп),Фм (Р1,Ум )}Т = 0 (39)
и, задав некоторое начальное приближение У® (по данным предыдущего шага по времени £ = £п) для последующих итераций в методе Ньютона из (36)-(38), имеем линейную относительно Усистему из 2М + 1 уравнений:
= 2(17),
3 = 0, 1, ■ ■ ■ , 3
В данном случае имеем
т т т
1 С1т-) ^ Р1т№$ч
(40)
dpi ’ dvim1 dpi
дф1 дф1 дф1
дрх ’ д^1т' др1
ТГ В'1т'У1т1 1тр1т1 1 1 =
др1т \
= {с1т ^1т^1т, ^1тр1т-) 1},
дух дщ д<Р1\_га „ а л А1
о •, о •, о г — \И1ш^1ш:И1шр1гп^] ^ др1т &01т др1 )
и непосредственной проверкой можно убедиться, что для нелинейной системы (36)-(38) её матрица Якоби дZ(Уj )/дУ?:
тхп =
щ3
Cii ±pii 0 0 .0 0 0
0 0 Ci2 ±Pi2 . .0 0 0
0 0 0 0 ci m ±PiM 0
(4 - km) -kipii 0 0 .0 0 -1
0 0 (- Ш -Ri2Pi2 . .0 0 -1
0 0 0 0 (cM - RiMviM) -Rmpm -1
вт kpn в 2v i 2 kPi2 . A mvi m ftiMPiM 0
является неособенной, поскольку
det(dZ(Yj )/dYj) = 0 .
В том случае, когда в ближайшем к узлу l конце xm = Xm какой-либо из входящих в этот узел ветвей m поток автомобилей является докритическим («сверхзвуковым»), то есть при vim > cim (Ai > 0,А2 > 0) искомые параметры Vim = {pim,vim} в таких окончаниях ветви m могут быть рассчитаны так же, как и в (36)-(38), однако условия совместности (24) следует заменить разностной аппроксимацией ударной адиабаты (19):
^im(vim) = (рПт (p?m) - рШ^^))*
х(1/р?г+1 - 1/plm) - (
vlm
п+Ц2 _ lm
)2 =0, (41)
так как при приближении к точкам ветвления обычно происходят многочисленные затормаживающие поток междуполосные перестроения и переход движения в закритический режим движения. Аналогичным образом следует поступать при
vim > Cim (Ai > 0,А2 > 0), Xm =0 для выходящей
из узла ветви, так как при интенсивном движении маловероятным является докритический режим движения непосредственно на выезде из перекрестка (во внутренних узлах графа). В этом случае использование (41) неявно предполагает допустимость разрывных волн ускорения («ударных волн разрежения»), что требует специальных исследований. Впрочем, возникновение таких ситуаций на перекрестках (узлах графа) с интенсивным движением маловероятно. Это замечание (об использовании (41) вместо условий совместности (24)) относится также к условиям замыкания граничных условий на входах в уличную сеть.
III. Численная реализация модели и результаты расчетов
Трудность реализации предлагаемой модели связана прежде всего со сложной структурой графа уличной сети, значительной неопределенностью определяющих параметров, констант «уравнения состояния» и процессов регулирования (многопараметричность задач), поэтому полномасштабная их реализация требует значительных вычислительных ресурсов и усилий по идентификации определяющих параметров («калибровка модели»). Структура искомого решения в таких системах может иметь очень сложный характер из-за взаимодействия многочисленных волн разгона — торможения друг с другом и с узлами графа, в том числе иметь разрывный характер. Это требует специальных исследований и предъявляет очень жесткие требования к выбору численных методов, пригодных для решения данного класса задач.
Составляющая основу данной модели система уравнений в частных производных (1), (2) име-
l
ет гиперболический тип и для численного решения таких уравнений существует достаточно много разнообразных по своим качествам численных методов. Однако отмеченные выше особенности потребовали специальных исследований и разработки новых монотонных схем высокого порядка аппроксимации ([17] и др.).
Описанная выше математическая модель была реализована в виде комплекса программ с использованием указанных численных методов и оттестирована на ряде задач. Эти тестовые расчеты показали высокую эффективность методов и работоспособность модели.
Определяющими функционирование описанной выше динамической модели уличного движения параметрами прежде всего являются характеристики графа уличной сети и уравнения состояния (2). Для оценки влияния неопределенности в их выборе (по имеющимся в литературе экспериментальным данным) были проведены тестовые расчеты автомобильного движения, примерно соответствующего условиям экспериментальных измерений [13, 14] с использованием уравнений состояния (13), (14) с константами из (29)--(32) (кривые а-й на рис. 2). Точное воспроизведение представленных в [13, 14] условий автодвижения принципиально невозможно из-за наличия на исследуемых участках пересечений с поперечными дорогами, по которым не приведено никаких данных, хотя их влияние на измеряемые параметры очевидно и может быть достаточно существенным, так как в местах пересечений обычно интенсивно происходят межполосные перестроения, что приводит к замедлению потоков, их уплотнению и даже к образованию заторов. Тем не менее по представленным в этих работах данным на входе (без учета участников движения с поперечных автодорог) были проведены расчеты и получены данные, качественно сопоставимые с экспериментальными данными на выходном измерительном пункте.
Рис. 9. Черными точками 1 показаны зависимости д(Ь,Хк) от р(Ь,Хк), полученные в различные моменты времени Ь. Темно-серыми точками 2 показаны те же зависимости для участка дороги, на котором реализуется свободный («сверхзвуковой») поток. Тонкие светлосерые линии 3 характеристики первого семейства, темно-серые 4 — ударные адиабаты
В одном из серии расчетов планарный дорожный граф состоял из двух участков, на входе первого участка задавался поток как возрастающая функция времени (33) (в соответствии с экспериментальным данными), а на выходе второго участка в качестве граничных условий задавались нулевые производные от плотности и скорости (34) (неотражающие условия, так как движение на втором участке соответствовало режиму «свободного потока). Во внутреннем узле графа одной из констант расчетной модели является коэффициент сопротивления Кт ^ 0 (36), который варьировался в довольно широких пределах: от нуля (когда такая граничная точка ничем не отличается от внутренних узлов разностной сетки) до очень больших его значений (соответствующих полному запиранию потока в этом месте).
На рис. 9 для одного из вариантов расчета в фазовой плоскости {д,р} с координатами Хк черными точками нанесены значения зависимостей д{Ь,Хк) от р(Ь,Хк), полученные в узлах разностной сетки первого участка дороги с координатами Хк в различные моменты времени £, то есть параметры, которые нанесены в виде экспериментальных точек на рис. 4. Темно-серыми точками (постепенно сливающимися в линию) нанесены аналогичные данные для второго участка, на котором реализуется свободный («сверхзвуковой») поток. Тонкими черными линиями нанесены те же, что и на рис. 4, характеристики первого семейства, а более толстыми светло-серыми линиями — ударные адиабаты для использовавшегося в расчетах уравнения состояния с константами из (30), поскольку сеточных узлов было 100, а моментов времени несколько тысяч, то различить можно лишь отдельные траектории в фазовой плоскости (самые толстые — черные линии). Однако четко просматриваются движения потока как близкие характеристическим направлениям (с «безударным» режимом торможения), так и движения с переходом от докритическо-го («сверхзвукового») режима к закритическому («дозвуковому») вдоль ударных адиабат. Для рассматриваемого здесь примера зависимость плотности автомобилей (черная кривая № 1), отнесенной к ее максимально допустимым значениям (пропускной способности первого участка), от координаты Хк, а также аналогичные распределения скорости (светло-серая кривая № 2) и потока (темно-серая кривая № 3) в один из моментов времени показаны на рис. 10. Эти два рисунка отчетливо демонстрируют, что закритический («синхронизированный») режим движения возникает в рассматриваемой здесь сугубо детерминистской модели в результате взаимодействия многочисленных волн разгона-торможения, возникающих и распространяющихся на участках дороги между локальными, замедляющими поток источниками (сужения дороги, пересечения с поперечными дорогами, различные замедляющие поток
маневры, неоднородность состава автотранспортных средств и стилей вождения и т.д.).
хъ 1
Рис. 10. Зависимости величин, отнесенные к максимальным значениям, от координаты хк: плотности автомобилей (черная кривая № 1); скорости (светлосерая кривая № 2); потока (темно-серая кривая № 3)
Рис. 11. Граф Садового кольца г. Москвы: 105 узлов; 212 ребер
Рис. 12. Результаты рассчитанных значений потоков [АТС / с] на Садовом кольце г. Москвы: слева для двустороннего движения; справа для одностороннего
Для проверки работоспособности разработанной вычислительной модели были сделаны численные расчеты поведения транспортной системы города Москвы внутри графа Садового кольца, имеющего 105 узлов и 212 ребер и представленного на рис. 11. В начальный момент времени
на входах в Садовое кольцо задавались постоянные значения потоков (33), полученные из данных экспериментальных замеров, сделанных в 2005 г. в центральной части г. Москвы, на выходах задавались неотражающие граничные условия (34). В результате расчетов были получены стационарные распределения потоков и скоростей движения автомобилей за 3-часовой интервал с момента начала движения. Затем граф данной транспортной системы перестраивался так, чтобы движение по Садовому кольцу становилось односторонним против часовой стрелки, и расчет проводился заново, стартуя с аналогичных исходных данных по значениям транспортных потоков на входах в транспортную сеть. В результате был получен эффект увеличения пропускной способности транспортной сети на 30%, что может быть вполне объяснимо, если исходить из того, что ширина дорог выросла в два раза при том, что средний пробег между любыми двумя точками увеличился меньше, чем вдвое из-за хорошей связности графа транспортной сети внутри Садового кольца. Результаты рассчитанных значений потоков [АТС / с] приведены на рис. 12 слева для двустороннего движения и справа для одностороннего.
Довольно подробный обзор по гидродинамическим моделям транспортных потоков имеется, например, в главах 2, 3 учебного пособия [18].
Работа поддержана грантами РФФИ № 10-07-00620-а, 08-01-00959-а, 08-07-00501-а,
08-07-00158-а, 10-07-00620-а, 11-01-00494-а, РГНФ № 08-02-00347, ПФИ ОМН РАН № 3, ПФИ Президиум РАН П-2. Работа проведена в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы.
Литература
1. Kholodov A.S., Kholodov Y.A. Computational models on graphs for the nonlinear hyperbolic system of equations // Proceedings of ASME 2004PVP Conference. — 2004. — V. 476. — P. 161-167.
2. Kholodov Y. A, Kholodov A.S., Kovshov N. V., Simakov S.S., Severov D.S., Bordonos A.K., Bapaev A.Z. Computational models on graphs for nonlinear hyperbolic and parabolic systems of equations // Proceedings of the III European Conference on Computational Mechanics. — 2006. — P. 2279.1-2279.19.
3. Морозов И.И., Холодов Я.А. Моделирование динамики транспортных потоков // Труды 51-й научной конференции МФТИ «Современные проблемы фундаментальных и прикладных наук». — 2008. — Т. 2. — С. 128-129.
4. Уизем Дж. Линейные и нелинейные волны. — М.: Мир, 1977.
5. Payne H.J. Models of freeway traffic and control. Mathematical Models of Public Systems // Simulation Council Proc. 28. — 1971. — V. 1. — P. 51-61.
6. Helbing D. Improved fluid: dynamic model for vehicular traffic // Phys. Rev. E. — 1995. — V. 51. — P. 3163-3169.
7. Aw A., Rascle M. Resurrection of «second order» models of traffic flow // SIAM Journal of Applied Mathematics. — 2000. — V. 60. — P. 916-938.
8. Greenberg J. M. Extensions and amplifications of a traffic model of Aw and Rascle // SIAM J. Appl. Math. — 2001. — V. 62, N. 3. — P. 729-745.
9. Zhang H.M. A non-equilibrium traffic model devoid of gas-like behavior // Transp. Res. B. — 2002. — V. 36. — P. 275--290.
10. Siebel F, Mauser W. Synchronized flow and wide moving jams from balanced vehicular traffic // e-print arXiv: physics /0509124 v2, 2006.
11. Kerner B.S. Introduction to modern traffic flow theory and control. The long road to three — phase traffic theory. — Springer, 2009.
12. Kerner B.S. and H. Rehborn. Experimental features and characteristics of traffic jams // Phys. Rev. E. — 1996. — V. 53, N. 2.
13. Kerner B.S. and H. Rehborn. Experimental properties of complexity in traffic flow // Phys. Rev. E. — 1996. —V. 53, N. 5.
14. Kerner B.S. and H. Rehborn. Experimental Properties of Phase Transitions in Traffic Flow // Phys. Rev. Let. — 1997. — V. 79, N. 20. -P. 4030--4033.
15. Kerner B.S. Experimental Features of Self-Organization in Traffic Flow // Phys. Rev. Let. — 1998. — V. 81, N. 20.
16. Годунов С.К., Роменский Е.И. Элементы механики сплошных сред и законы сохранения. — М.: Наука, 1998.
17. Холодов А.С., Холодов Я.А. О критериях монотонности разностных схем для уравнений гиперболического типа. // Журнал выч. математики и мат. физики. — 2006. — Т. 46, № 9. — С. 1560-1588.
18. Гасников А.В., Кленов С.Л., Нурмин-ский Е.А., Холодов Я.А., Шамрай Н.Б. Введение в математическое моделирование транспортных потоков; учебное пособие. С приложениями Бланка М.Л., Гасниковой Е.В., Замятина А.А. и Малышева В.А., Колесникова А.В., Райгородско-го А.М. — М.: МФТИ, 2010.
Поступила в редакцию 15.10.2010.
