CC BY
30
УДК 519.214.6, 519.722, 519.876
А.В. Гасников, Е.В. Гасникова
Московский физико-технический институт (государственный университет)
О возможной динамике в модели расчета матрицы корреспонденций
(А.Дж. Вильсона)
Приводится основной аппарат, необходимый для исследования динамики макросистем при больших значениях времени. В основе динамики лежит эргодическая марковская цепь с огромным числом состояний. При больших значениях времени распределение макросистемы по макросостояниям будет близко к стационарному. С ростом размерности макросистемы (количества состояний марковской цепи) стационарное распределение будет концентрироваться в окрестности наиболее вероятного макросостояния, которое и принимается за равновесное для данной макросистемы. В качестве примера применения описанного формализма приводится вывод статической гравитационной модели расчета матрицы корреспонденций (одной из наиболее популярных в приложениях), исходя из «разумной» (индивидуально выгодной) динамики обменов местами жительства.
Ключевые слова: эргодическая теорема, функция Ляпунова, энтропия, гравитационная модель расчета матрицы корреспонденций, концентрация инвариантной (стационарной) меры, канонический скейлинг, динамика квазисредних, условие динамического равновесия, принцип детального равновесия.
I. Введение
После работ ряда крупных биологов первой половины прошлого века, а также работ Е.Т. Джейнса (конца 50-х годов XX века) [1], А.Дж. Вильсона (конца 60-х годов ХХ века) [2], И. Пригожина с коллегами, Г. Хакена (70-е годы XX века) [3, 4] в литературе достаточно прочно укрепилась концепция о плодотворности перенесения термодинамического формализма (см., например, [5-14] и цитированную там литературу) на различные макросистемы (в частности, встречающиеся в экономике, биологии, социальной сфере [2-4; 15-24]). В России систематические исследования в этом направлении были предприняты Л.Н. Розоноэром в начале 1970-х [25] (см. также [26-34] и цитированную там литературу). Упомянутая концепция часто используется для нахождения равновесия макросистемы. А именно, по аналогии с феноменологической термодинамикой вводится вероятностное распределение на множестве состояний, в которых может пребывать макросистема. Такое распределение может, например, совпадать с инвариантной мерой эргодической динамической системы, порождающей рассматриваемую макросистему [11], или с финальным (равным стацио-
нарному) распределением эргодического (например, марковского) случайного процесса, порождающего рассматриваемую макросистему [35-40]. Если размерность макросистемы увеличивается, то, как правило, распределение сосредотачивается в окрестности наиболее вероятного макросостояния7. Таким образом, с ростом времени наблюдения за макросистемой и размерности макросистемы следует ожидать нахождения макросистемы с большой вероятностью в малой окрестности наиболее вероятного макросостояния вне зависимости от того, в каком состоянии макросистема находилась сначала (иначе говоря, большую часть времени (иногда, и просто, на больших временах) макросистема будет пребывать в малой окрестности наиболее вероятного макросостояния). Естественно поэтому под равновесием макросистемы понимать наиболее вероятное макросостояние. Задача нахождения наиболее вероятного макросостояния часто сводится (асимптотически по размерности системы) к задаче максимизации энтропийно подобного функционала при ограничениях (в термодинамике таким образом можно получить статистики Больцмана, Ферми-Дирака, Бозе-Эйнштейна [1, 5]). Подробнее о приложениях этой концепции см., например, [1, 2, 30, 33; 41-45]. 8
7Заметим, что отмеченное обстоятельство (концентрация) может быть по-разному обосновано (как правило, достаточно элементарных комбинаторных соображений и формулы Стирлинга [1, 2, 30, 33].
8Укажем некоторые часто встречающиеся в приложениях [41-45] формализмы, также приводящие к задачам оптимизации энтропийно подобных функционалов: принцип наибольшего правдоподобия (при оценке неизвестных параметров по простой выборке); принцип максимума апостериорной вероятности; наименьшее отклонение в смысле расстояния Куль-бака—Лейблера (энтропийное расстояние); принцип наименьшей неопределенности (энтропия — мера неопределенности) в теории информации (рассуждения опираются в ряде случаев на теорему Шеннона—МакМиллана—Бреймана). Важно также отметить, что энтропийно подобный функционал часто является функцией Ляпунова для динамической системы (например, системы обыкновенных дифференциальных уравнений, итерационного процесса, уравнения (-ий) в частных производных эволюционного типа и т.п.), порождающей рассматриваемую макросистему [12; 17—20]. Пожалуй, наиболее ярким примером этого тезиса является кинетическая теория (Л. Больцман [8]).
II. Возможная динамика, приводящая в асимптотике (по времени) к статической модели А.Дж. Вильсона расчёта матрицы корреспонденций
Рассмотрим для начала более простой пример, иллюстрирующий формализм, описанный во введении.
Пример 1 (кинетика социального неравенства [23, 46]). В городе живет N ^ 1 (например, 10 000) пронумерованных жителей. У каждого *-го жителя есть в начальный (нулевой) момент времени целое (неотрицательное) количество рублей в*(0) (монетками, достоинством в один рубль). Со временем пронумерованные жители (количество которых не изменяется, так же как и суммарное количество рублей) случайно разыгрывают свое имущество. Пусть в момент времени £ ^ 0 г-й житель имеет к рублей, а 1-й житель — т рублей. Тогда Рк-т(£)Д£ + о(Д£) есть вероятность того, что жители с номерами г и I (1 ^ г < I ^ N) встретятся и попробуют разыграть один рубль по следующему правилу: с вероятностью 1/2 житель с большим номером отдаёт 1 рубль (если, конечно, он не банкрот) жителю с меньшим номером и с вероятностью 1/2 наоборот. Будем считать, что Р^,т (£) = кN-1 (к > 0). При этом «в среднем» в единицу времени осуществляется кN/2 встреч, т. е., скажем, при к = 1 в единицу времени каждый житель с вероятностью, большей 1/2, с кем-то должен встретиться. Приблизительно такую постановку задачи в конце XVIII века предложил известный итальянский экономист Вильфредо Парето, чтобы объяснить социальное неравенство.
Пусть с5(Ь) — доля жителей города, имеющих ровно в рублей в момент времени £ (заметим, что с5(£) — случайная величина). Пусть
Тогда по эргодической теореме для конечных однородных марковских цепей (см. [18, 19; 35-38] и упражнение ниже): 9
V д = 0, Я 3 Хд > 0, Тд = О^) : V £ > Тд ^
Р
1
< в
0, ..., д > 0,999
Се-5/5
где С определяется из условия нормировки:
ЕСе~
5 = 0
= 1, то есть С ж
в
Скорость сходимости оценивается сверху, исходя из оценок в доказательстве эргодической теоремы для однородных марковских цепей с конечным числом состояний.10 Как показывают численные эксперименты, оценка О^) точная. Так, если в городе 10 000 жителей и единица времени — день, то при начальном «социальном равенстве» с вероятностью, близкой к единице, через 20-30 лет (при к = 1) установится «социальное неравенство». Заметим, что описанный выше случайный процесс обратим во времени. Однако наблюдается необратимая динамика с3(£). Но в таком случае можно удивляться также и тому, что газ, собранный в начальный момент в одной половине сосуда, с течением времени равномерно распределится по сосуду.
Приведем отчасти схожую постановку задачи (та же самая мера будет концентрироваться), восходящую к В.П. Маслову [33]. Ниже приведен фрагмент его интервью 2009 года, посвященного объяснению финансового кризиса 2008 г.
В.П. Маслов: — Поясню на знаменитом трюке Коровьева-Фагота. Помните булгаковского героя, который разбрасывал в варьете червонцы? Понятно, что кому-то досталось больше купюр, кому-то меньше, а кто-то вообще остался ни с чем. Вопрос: если купюр миллион, то сколько должно быть зрителей, чтобы ни один не остался без банкноты? Вроде очень неопределенная задача, не имеющая однозначного решения. И тем не менее ответ
9Эргодическая теорема используется для нахождения распределения случайных величин Се() на больших временах. Далее используются законы больших чисел или, другими словами, явление концентрации инвариантной (стационарной) меры, о котором мы подробнее поговорим в следующем примере. Точнее не само это явление, а его следствие о том, что «хорошие» (например, липшицевы) функции на «хороших» компактных пространствах с мерой большого числа измерений почти везде близки к медиане [47].
10 Необходимо (для простоты формулировок, в рамках этой сноски считаем время дискретным) асимптотически (по размеру макросистемы) оценить второе по величине модуля собственное значение матрицы переходных вероятностей — инфини-тезимальной матрицы (первое собственное значение, которое для неотрицательных матриц часто называют числом Фробе-ниуса—Перрона, равно единице, поскольку матрица стохастическая — все элементы неотрицательны и сумма всех элементов в любой строке равна единице), определяющее основание геометрической прогрессии, мажорируемой последовательность норм отклонений текущего состояния от стационарного в различные моменты времени [35-38]. Здесь нельзя не упомянуть о том, что в этом направлении за последние несколько десятков лет произошла определенная революция [48], которую можно пояснить рассмотренным примером 1. Несложно проверить, что число (макро-) состояний марковской цепи в этом примере растет быстрее, чем экспоненциально с ростом N. В то время как по прошествии всего лишь 0(^) тактов распределение цепи будет уже довольно близко к финальному (предельному) = стационарному (инвариантному). Таким образом, если возникает потребность быстро сгенерировать дискретные случайные величины, которые могут принимать огромное число значений, то в ряде случаев удается подобрать такую марковскую цепь, которая быстро «выйдет» на стационарный режим, соответствующий желаемому распределению. Несколько интересных примеров в этом направлении (модель Изинга и др.) собрано, например, в современном курсе марковских случайных процессов [38]. Заметим, что при оценках второго по величине модуля собственного значения активно используется уже упоминавшийся принцип концентрации меры (см. [48, 49] и цитированную там литературу, а также статью А.В. Колесникова).
N
Я
i=1
я
м'
с
есть: примерно квадратный корень из миллиона, то есть тысяча зрителей.
Точнее говоря, как следует из вышенаписанно-го, с вероятностью, близкой к 1, доля банкротов будет равна примерно 1/е ~ 0,001, поскольку по условию вЖ = Б = 106 и N ~ д/б*. Поэтому количество банкротов с вероятностью, близкой к 1, незначительно отличается от N/8 ~ 1.
Пример 2 (модель расчета матрицы корреспонденций [2]). В некотором городе имеется п районов, ^ > 0 — число жителей г-го района, Шу > 0 — число работающих в ]-м районе (число рабочих мест), ху(£) ^ 0 — число жителей, живущих в г-м районе и работающих в ] -м в момент времени £ ^ 0. Со временем пронумерованные жители (количество которых не меняется11 и равно N = ^П=1 ^ = ^П=1 Шу) меняют места жительства (квартиры). Считается, что отмеченные изменения могут происходить только за счёт обмена квартирами, то есть
хі3 (^) ^ 0- ^ ' !'/у (^) — ^і-
3=1
^Хіз (І) — - І- І = I-
(А)
Пусть в момент времени £ ^ 0 г-й житель живет в к-м районе и работает в т-м, а 1-й житель живет в р-м районе и работает в д-м. Тогда Рк т-р д(£)Д£+о(Д£) есть вероятность того, что жители с номерами г и I (1 ^ г < I ^ N) поменяются
квартирами в промежутке времени (£,£+Д£). Естественно считать, что вероятность в единицу времени обмена местами жительства зависит только от мест проживания и работы обменивающихся.
Например, можно считать, что время, потраченное в пути от района г до района ], есть Ьу ^ 0 (вместо в приводимую ниже формулу также осмысленно подставлять ^ 0 — расстояние от района г до района j), а
Рк,
7(І) —
= Р N ехр(7 • ^т + £рд - (£рт + tkq))/2) > 0,
где рк > 0, ^ > 0. Тогда по эргодической теореме для конечных однородных марковских цепей (см. [18, 19; 35-38]):
V {ху }”Ду=1 е (А) -
— Иш р(ху (г) = ху, г, j = 1,..., п) =
-і
ехр(-7^3Ху) • (хіз!) 1 ^
=1,3=1
йєі
і3
-їП,П \
}і=1,3 = 1)-
где статсумма ^ находится из условия нормировки получившейся пуассоновской вероятностной меры [50]. Распределение р({ху}П==1 у =1) на множестве (А) сконцентрировано при N ^ 1 (см. ниже) в окрестности наиболее вероятного значения
11 Поскольку мы будем следить за системой на больших временах, то сделанное предположение кажется неестественным. Заметим, однако, что если «номер жителя» передаётся его потомкам (номер папы передаётся сыну, номер мамы — дочери), то предположение о постоянстве состава жителей выглядит разумным в первом приближении. Здесь мы также пренебрегаем миграционными потоками (город изолирован).
12Нетрудно заметить, что задача поиска наиболее вероятного состояния асимптотически (по п) эквивалентна задаче максимизации энтропийного функционала (воспользовались формулой Стирлинга):
1п р({хіз ^.^
1) Хіі 1п(х^' /е) - 7 ^2 їіз Хіі + ООПв^п
г=1,^=1 і=1^=1
на множестве (А). Поскольку функционал строго вогнутый и решение задачи максимизации без предположения целочис-ленности {xij}= ^_1, считаем таким, что х^ ^ 1 (так как п ^ 1), то ограничение «целочисленности» является для данной задачи несущественным и применение асимптотической формулы Стирлинга было законным. Обратим внимание также на то, что задача максимизации энтропийного функционала на множестве (А), то есть по принципу Лагранжа [51] (балансовые ограничения (А) перенесли в функционал).
L({xij}п_1^=1; ЛЬ,ЛШ) = - ^2 Х^ 1п(хі,/е) - 7 ^2 tijХ^ +
і = 1^' = 1 і = 1^' = 1
+£ л%
і=1
І п \ п / п \
- ЬА + £ ЛТ • ( £ Хij - ^ І
\j-_1 ) j=l \і=1 )
8Ир
0 ,i,j = 1
имеет, и притом единственное решение {Х^}
п,п 1^'
дL({xij }
п,п 1^'=
которое определяется из системы: ____ 1- . I \Ь
; ЛЬ,ЛУ )/дх^ = - 1п Хij + Х% + ЛУ - = 0, = 1,...,п.
Приятной особенностью класса задач энтропийно линейного программирования (задач максимизации энтропийно подобных функционалов, при наличии линейных ограничений в виде равенств и неравенств на неотрицательном ортанте) является явная (легко выписываемая) зависимость решения прямой задачи, через двойственные переменные. Поскольку количество ограничений (количество двойственных переменных — множителей Лагранжа), как правило, намного порядков меньше числа прямых переменных, то эффективные численные методы отыскания решений базируются на решении двойственной задачи, представляющей собой задачу минимизации выпуклой функции на неотрицательном ортанте [30, 45, 52]. Можно показать, см. [52], что многие популярные в литературе [30, 45] численные методы решения этой двойственной задачи являются барьерно-мультипликативными аналогами квазиградиентных итерационных методов. В частности, к ним относится один из первых методов: «метод балансировки», заключающийся в подстановке прямых переменных как функций двойственных, в
п. п
п. п
п.п
п.п
п
г * -1 п,п
{ху ^=1 у=1, которое находится, как решение задачи энтропийно линейного программирования:
п,п i=1,j = 1
ху 1п(ху/е) + 7
Е
х1.
Перейдем теперь к исследованию полученного стационарного распределения вероятностей р{{ху }П=П у=1) на макросостояниях {ху}п=п1 ,у=1 е (А). Если13
1,у=1
N — п“+1, и,Шу —
па (а > 1), п ^ 1,
(1)
{х^3 }г= 1,^ = 1
Решение задачи (1) можно представить как
ху = ехр( —1 — хк — у — ),
где множители Лагранжа (двойственные переменные ) {Хк}”=1
и {Х^ }п=1 определяются из равенств системы (А). На практике мы имеем информацию о {Ь,,,Ш,,}п=1 и {1у}п=п1 у=1. Решив задачу (А), мы найдем
{хт({ЬШ}п=1; {Цу}”=п1 ,у=1)}п=п1 ,т=1.
Такой способ расчета матрицы корреспонденций в литературе часто называют гравитационной моделью:
xij АiBjLi^Vjехр( ),
где {Аi}п=1 и {Ву}п=1 определяются из соотноше-
V г, j = 1,..., п — иу = т > 0, то распределение вероятностей р({ху }п=п у=1) на множестве (А) сконцентрировано в О(п(а—1)/2) окрестности (почему именно в такой окрестности, будет показано ниже) наиболее вероятного значения:
х*у « ЬШу /N — п'
а—1
1.
(2)
ний [2, 30, 32]:
А= \^2 В у Шу ехр( )
у=1
1
1
которое ищется с помощью метода множителей Лагранжа [2, 51]. Формулу (2) можно интерпретировать как наличие у районов «потенциалов притяжения» [2]:
Ьг/^[Ы ~ ехр(Л^), г = 1, п
ТУ^/у/Ы - ехр(Л^), ] = 1, п,
произведение которых даёт пассажиропоток х* у, г, 3 = 1, ..., п (асимптотически совпадающий с математическим ожиданием и медианой).
Исследуем теперь, следуя [1, 2, 30], явление концентрации стационарного распределения р(.{хну }п=" у=1) в окрестности наиболее вероятного
п п
метим, что из определения {ху }i=l у = 1 (см. также
Ву
Е АЛ ехр(-^ку)
^=1
значения {х* у }п=п у=1. Для этого прежде всего заметим, что из опре сноску 12) следует
V {ху}п=п1 ,у=1 е (А)
п п
а 1п р({ху }п=п1, 7=1)
i=1 , у=1
°цп=1,3 = 1/
х (xij — х\у) ^ 0.
V{xгз }п=п1 у=1 е (А) 3 в е [0,1]:
1пр({ху }п=п1 ,,=1) < 1пр({х*3 }п=п1 ,,=1) +
^ д2 \пр{{х*ув + хц • (1 - 0)}^" =1)
+ 2^ ------------------------------~-----------------—х
*=1,у = 1
дх2
дху
{ХЮ - Х*з)2
Описанная модель (точнее говоря, рассчитанная по этой модели матрица корреспонденций) активно используется в теоретико-игровых моделях (например, базирующихся на принципах Дж.Г. Вар- Поэтому
дропа) равновесного распределения потоков [32]
(см. также главу 1 в пособии [54]). Подробнее об этих моделях речь пойдет в статьях Е.А. Нурмин-ского и Н.Б. Шамрай, В.И. Швецова, Е.В. Гасни-ковой и Ю.В. Дорна из этого же сборника. Один из возможных способов определения времени в пути, в зависимости от загрузки ребра, предложен в статье М.Л. Бланка. Иной подход к выводу гравитационной модели имеется в статье П.П. Бобрика.
Заметим также, что задача (1) по своим свойствам Но
очень похожа на транспортную задачу (см. статью А.В. Колесникова).
систему ограничений (в методе предполагается, что есть ограничения только в виде равенств), и разрешение полученной системы (размерность которой как раз равна числу двойственных переменных) относительно двойственных переменных. Для рассматриваемого далее частного случая V г, ^ = 1, ..., п = т > 0 все это можно сделать аналитически и в результате
получить формулу (2). Отметим здесь также эффективность сепарабельных (функционал декомпозируется в аддитивную сумму функций одного аргумента) алгоритмов типа Нестерова—Немировского для задач энтропийного программирования, возникающих при нахождении равновесий макросистем [53].
13Отметим, что хотя в этом случае динамика рассматриваемой нами макросистемы обратима по времени (так же, как и в примере 1), макросистема (в каком бы состоянии она не находилась в нулевой момент времени) по прошествии достаточно большого времени окажется в малой окрестности равновесного макросостояния (характеризующегося наибольшим из возможных значением энтропии) и будет в дальнейшем пребывать в этой окрестности подавляющую часть времени. Схожая ситуация имеет место и в статистической физике (см., например, [8, 11, 14, 18, 19, 27]).
а2
дх2-
дху
(1п р({хго }п=п1,у=1))
и
д2
дх2 dxij
n,n
n,n
= l,j=
xij ln xi
Следовательно, приходим к «неравенству о концентрации меры»:
V M > 0, V {xij }
n,n = 1,j =
Є (А):
n,n
n,n
=l,j=l
(Xjj - X*j )2
2тах{жу,.г-*.}
'P({xij }
) < e-Mp({x*j}=1j=l).
Из этого неравенства имеем результат о концентрации распределения р({х,,у}п=п1 у=1) на множестве (А) в О(п(а—1)/2) окрестности наиболее ве-
п п
роятного значения {х}i=l у=1:
З A> 0 : lim P (V i, j = 1,
Замечание (о других возможных подходах к исследованию концентрации стационарного распределения). Один из способов восходит к методу вычисления математических ожиданий Дарвина-Фаулера [2, 5, 9] (метод производящих функций и анализ их асимптотического поведения методом перевала) — в этом случае концентрация наблюдается в окрестности математического ожидания (интересные приложения этого метода в комбинаторике имеются, например, в [55], см. также статью А.А. Замятина и В.А.
Малышева в этом сборнике и конец следующего пункта). Исследование концентрации в окрестности математического ожидания можно также приводить, например, используя предельные меры [56], метод канонических ансамблей [57] или обобщённую схему размещения [58], нашедшие применения к задачам асимптотической перечислительной комбинаторики,14 к исследованию случайных матриц и уравнений, к изучению статистических свойств группы перестановок с приложениями к теории разбиений (диаграммам Юнга) и асимптотической теории чисел, а также к теории предельных форм. К методу производящих функций также тесно примыкают метод моментов [58], метод пуассоновской и гауссовской аппроксимации (метод локальной предельной теоремы) [7, 58].
Другой способ восходит к принципу концентрации А. Пуанкаре и П. Леви, получившему дальнейшее развитие в работах В.Д. Мильмана и др. [59] — в этом случае концентрация наблюдается в окрестности медианы.15 В заключение краткого обзора
14В частности, к теории случайных графов (компьютерные, транспортные сети — вопросы надёжности и т.п. [24, 39], см. также статью А.М. Райгородского в этом сборнике).
15 Этот принцип также нашёл широкие применения в асимптотической перечислительной комбинаторике; в качестве достаточно известного примера можно упомянуть неравенство М. Талаграна и его приложения к изучению макросвойств (связность и т.п.) случайных графов [60] (см. также статьи А.М. Райгородского и А.В. Колесникова в этом сборнике).
методов исследования концентрации меры упомянем теоремы тауберова типа [61] и мартингаль-ные неравенства [60]. Ввиду всего вышесказанного важно отметить, что поиск наиболее вероятного распределения — это особенность работ, в которых изучаются равновесия макросистем. Но также важно заметить, что наиболее вероятное распределение в содержательных задачах асимптотически (по размеру системы) эквивалентно математическому ожиданию (в работе [2] соответствующие выкладки проделаны на примере модели расчета матрицы корреспонденций А.Дж. Вильсона, см. также следующий пункт) и медиане [59].
Замечание (модели Д. Бернулли—Лапласа, П. и Т. Эренфестов). Важно обратить внимание на то, что описанный способ изучения равновесных состояний макросистем применим к достаточно широкому классу макросистем (модель Д. Бернулли-Лапласа, модель Эренфестов, круговая модель М. Каца и др. (см., например, [8, 27, 62] и цитированную там литературу)), в том числе встречающихся в экономике, биологии, социальной сфере [3, 4; 15-24].
III. Общая схема исследования равновесий макросистем
Ниже приводится (во многом под влиянием работ [18, 19]) общая схема, в которую «ложатся» примеры 1 и 2.
Предположим, что некоторая макросистема может находиться в различных состояниях, характеризуемых вектором п с неотрицательными целочисленными компонентами (скажем, в модели «кинетика социального неравенства» щ(1) — количество жителей города, имеющих в момент времени £ ^ 0 г рублей). Будем считать, что в системе происходят случайные превращения (химические реакции). Пусть п — п — а + 3, (3,3) е I — все возможные типы реакций, при этом (3,3) е I ^ (3,3) е I. Введем, следуя Березину-Клесову (1978) [19], интенсивность реакции:
A(a , в)(п) = Л
M
где
Kа
[ щ • ... • (щ — аi + 1),
i:ai>0
^ 0 — константы реакции (в химической кинетике — постоянные, а в социодинамике (В. Вайдлих, 2000 [20]) — необязательно); при этом часто считают 1 п = М, т. е. Х(йД)(п) — вероятность осуществления в единицу времени перехода п — п—3+3: в единицу времени равновероятно выбираются («приближение среднего поля») любые два жителя города и в зависимости от того,
І
x
ij
n,n
n
x
в каких состояниях они находились, «случайно» переводятся (разыгрывают один рубль) в новые состояния. На макроуровне все это соответствует принципам химической кинетики (закон действующих масс Гульдберга-Вааге, 1864 [18]).
Предположим теперь, что множество .1 не зависит от М, и в начальный момент времени для любого І существует предел
сі(0) = Ііт пі(0)/М.
М
Тогда (Малышев-Пирогов-Рыбко, 2004 [19, 63]) в произвольный момент времени Ь > 0 и для любого і существует предел по вероятности (заметим, что пі(Ь) — случайные величины, тем не менее сі(Ь) — уже не случайные величины):
Несложно проверить, что 3, удовлетворяющее (4), является неподвижной точкой системы (3). Более того, если существует хотя бы одно положительное решение системы (4), то тогда все неподвижные точки системы (3) удовлетворяют условию (4), также называемому условием унитарности, или условием Штюкельберга-Батищевой-Пирого-ва. Покажем, во многом следуя В.В. Веденяпину [18], что в этом случае (существования хотя бы одного положительного решения (4)) траектория (3) сходится к неподвижной точке (какой именно, зависит, вообще говоря, от «точки старта»).17 Для этого, следуя второму методу Ляпунова, введём (минус) энтропию:
Н
^2Сі • (Іп(с/і) - 1)
сі(Ь) =' Ііт пі(Ь)/М.
М
Описанный выше приём называется каноническим скейлингом. В результате такого скейлинга приходим к «динамике квазисредних» (терминология В. Вайдлиха [20]):
ё,с.
(3)
(аф)Е.]
Эти же уравнения можно получить и по-другому. А именно, как приближенную динамику средних с(£) = Е[пД^/М]. Приближенную в том смысле, что при выводе (1) используется приближение
Г(8(I)) « Е[Г(т(1)/М)]
для «достаточно хороших» функций Г (например, полиномов).
Это верно в случае пикообразного распределения п^).16 Пусть существует хотя бы одно положительное решение 3 системы (условие динамического равновесия, В.В. Веденяпин (2001) [18], Малышев-Пирогов-Рыбко (2004) [19, 63], Бати-щева-Веденяпин [64]):
Е
0:(а,0)eJ 3
_Е
p:(а,p)ЄJ
4вз
3-
(4)
где К|(с) = Ка, а а — произвольный, но такой,
что 3 3: (а.,/3) е I.
Частным, но часто встречающимся в приложениях, случаем условия (4) является условие детального равновесия:
и покажем, что она является функцией Ляпунова для системы (3). Посчитаем полную производную Н в силу системы (3):
% = Е
(а,р)еЗ 3
ІпП у
ві — аі
- (в - аі)
+
р 11 ^3 *3
(аф)ЄІ 3
Е К!Ш у? •ув
уа
И^' у3
(аф)Е-1 3
ві — аі
где введено обозначение уі = сі/(і. Заметим, что
Е К!Ш’у? = Е К'Ш(в'у? =
(а,/3)€3
(а,Р)€3
Е К1П (
(аф)Е-І 3
Таким образом,
п
—в'
3
а' в'
? у в х
(а,Р)€3
1пП у? вг
Цу?—в' + 1
К|П 3 = КрП (в'.
поскольку и 1п и — и + 1 ^ 0 при и > 0, и равенство достигается в одной точке и =1.
Оказывается (Малышев-Пирогов-Рыбко [19, 63]), что условие (4) можно проинтерпретировать
16Заметим, что этот переход и возможность его использования в правой части равенства (3) нуждаются в строгом обосновании (и далеко не всегда правомочны). В качестве примера укажем популярный в литературе [3; 16—22] марковский процесс «рождения-гибели» (приводящий к системе уравнений «хищник-жертва»), для которого «флуктуации играют решающую роль, качественно меняя выводы макроскопического анализа».
17Стоит заметить, что аттрактор системы (3) с постоянными коэффициентами реакции, по-видимому, в общем случае может быть сколь угодно сложным множеством [19, 65].
X
X
как условие инвариантности пуассоновской меры [50]:18
М(п) = П ^ —/п^:
%
где Хi = £* М, а 3* — произвольное решение (4); относительно предложенной стохастической марковской динамики. Эта мера экспоненциально быстро концентрируется, с ростом щ, в окрестности наиболее вероятного состояния (также удовлетворяющего условию (4)), которое и принимается за положение равновесия макросистемы. Задача поиска наиболее вероятного макросостояния асимптотически эквивалентна задаче максимизации энтропийного функционала (воспользовались п\ = у/2тг(п/е)п(1 + о(1)) — формулой Стирлинга):
Тогда по формуле Коши:
11
Е [пР1 •
■пр2]
П
— (1к — 1
Z (2пг)т 1
ф^1 • ... • с!гтх д
1п ^1 дА1г
Рч
Z =
1
(2пг)т
()й^1 • ... • ё,гт-
П
Г (3; А)
Лк — 1 I • Г(3; А)
Е
Здесь математическое ожидание Е[пР1 •...•пр^] считается по вероятностной мере, порожденной мерой Пуассона ^(п), а интегралы по dzk берутся в комплексной плоскости по замкнутым контурам, охватывающим точку ноль. Используя метод перевала [67], асимптотически оценим математическое ожидание:
Е [пР1
■пр%]
1
на множестве, заданном (как правило, линейными) ограничениями — законами сохранения (интегралами движения для (3)).
Действительно, будем считать, что ограничения (законы сохранения) задаются СЛАУ Ап = d, где А = ||Ам|| — матрица максимального ранга (к = 1, ..., т). Обозначим через А множество неотрицательных целочисленных векторов п, удовлетворяющих Ап = d. Тогда равновесие п* находится как решение задачи Е(п) — шахп^л (поскольку функционал строго вогнутый и считаем п* ^ 1, то целочисленностью переменных можно пренебречь). Используя принцип Лагранжа, можно показать, что решение этой задачи представляется в виде
щ(у*) = Хi ехр
п
Г(3*; А) д ^
(1п z1)
Рч-
дА
1га
Г (3*; А)
где «точка перевала» 3 определяется как решение системы:
Zk • (дГ(3; А)/дzk) ~ dk • Г(3; А), к =1, ..., т.
В частности,19
Е[п] « Х П(zk А, D[ni] « Х П(zk)Аы,
k k
где 3* определяется как решение системы уравнений:
ЕА^ • <
,Аы
где двойственные переменные (множители Лагранжа) у* определяются из системы уравнений Ап (у) = 3.
Приведем, во многом следуя [2], другой путь (восходящий к Дарвину-Фаулеру), по которому можно прийти к аналогичным формулам.
Для этого введем производящую функцию:
Г(3;А) = Е м(п)П
>0
Очевидна связь «точки перевала» 3* с двойственными переменными у*: z*к = exp(yk).
IV. Заключение
В статье обсуждается концепция равновесия макросистемы. Приводятся различные подходы к обоснованию следующего принципа: равновесие = наиболее вероятное макросостояние инвариантной (стационарной) меры динамической системы (марковского процесса), порождающей исследуемую макросистему. Рассматриваются примеры конкретных макросистем. В частности, один из примеров «объясняет» популярную в приложениях модель А.Дж. Вильсона расчета матрицы корреспонденций.
Повторим в заключение описанную в статье схему.
18 Пуассоновская мера также часто встречается (при наиболее естественных условиях) в качестве притягивающей инвариантной меры, если осуществляется термодинамический предельный переход (число агентов и состояний стремятся, сохраняя пропорции, к бесконечности) [66].
19Обратим внимание на то, что получилось: Е[п] ~ ^[п^] ^ 1 — это означает концентрацию распределения случайной величины щ в \/ 1?[пг] окрестности своего математического ожидания п-; \.
■,-Хг
ПЕ
i п^ > 0
1
X
X
k
1. Макросистема состоит из огромного числа пронумерованных агентов, каждый из которых может находиться в одном из возможных состояний. Число состояний, как минимум, на несколько порядков меньше числа агентов (иногда можно и без этого требования). Распределение агентов (с учетом их номеров) по состояниям будем называть микросостоянием, а без учета номеров — макросостоянием.
2. Задана марковская динамика распределения агентов по состояниям, в основу которой на микроуровне положена равноправность агентов одного типа (в приближении среднего поля) и заранее прописанные возможности случайных превращений (переходов) агентов (химические реакции): равновероятно выбирается агент и в зависимости от того, в каком состоянии он находится, «случайно» переводится в новое состояние. Аналогично рассматриваются парные взаимодействия и взаимодействия, в которых участвует большее число агентов. На макроуровне это соответствует принципам химической кинетики (Гульдберг-Вааг, 1864).
Предполагается, что из любого возможного макросостояния можно перейти согласно такой динамике в любое другое (характерное время такого перехода определяет скорость сходимости к равновесию). Также считается, что описанная динамика имеет макрозаконы сохранения. Соотношения (как правило, линейные) между макровеличинами, которые не меняются со временем.
Пусть выполняется условие: динамика задана линейной полугруппой (однородность), динамика «обратима» (детальный баланс, условие динамического равновесия).
Тогда эргодическая марковская динамика приводит на больших временах к стационарной (инвариантной) пуассоновской (сложной) мере (прямое произведение распределений Пуассона) на пространстве макросостояний. Эта мера экспоненциально быстро концентрируется, с ростом числа агентов, в окрестности наиболее вероятного макросостояния, которое и принимается за положение равновесия макросистемы. Задача поиска наиболее вероятного макросостояния асимптотически (по числу агентов) эквивалентна задаче максимизации энтропийного функционала (воспользовались формулой Стирлинга) на множестве (как правило, аффинной структуры), заданном ограничениями — законами сохранения. Этот же энтропийный функционал (взятый со знаком минус) возникает как функция Ляпунова динамики, полученной в результате канонического скейлинга исходной марковской динамики. Отыскание предельной неподвижной точки (этой динамики), в которую придет система, сводится к решению той же самой задачи энтропийно линейного программирования. Приятной особенностью такого класса задач является явная (легко выписываемая) зависимость решения прямой задачи через двой-
ственные переменные. Поскольку число ограничений, как правило, на много порядков меньше числа прямых переменных, то эффективные численные методы базируются на решении двойственной задачи: задачи минимизации выпуклой функции.
Авторы выражают благодарность за ряд ценных замечаний В.В. Веденяпину, И.А. Лубашев-скому, В.А. Малышеву, Е.А. Нурминскому, С.А. Пирогову, Ю.С. Попкову, И.Г. Поспелову, В.Н. Разжевайкину, С.В. Чуканову, Н.Б. Шамрай, А.А. Шананину, В.И. Швецову.
Работа поддержана грантами РФФИ № 08-01-00959-а, 08-07-00501-а, 08-07-00158-а, 10-01-00321-а, 10-07-00620-а, 11-01-00494-а, РГНФ № 08-02-00347, ПФИ ОМН РАН № 3, ПФИ Президиум РАН П-2. Работа проведена в рамках реализации ФЦП «Научные и научно-педагогические кадры инновационной России» на 2009-2013 годы (мероприятие 1.2.1, НК-15П, П949; мероприятие 1.3.1, НК-215П, П1490).
Настоящая статья представляет собой запись нескольких лекций, прочитанных первым автором студентам МФТИ в весеннем семестре 2009/2010 учебного года в рамках курса по выбору «Математическое моделирование транспортных потоков».
Литература
1. Jaynes E.T. Probability theory. The logic of science. Cambridge: Cambridge University Press, 2003; Papers on probability, statistics and statistical physics. Dordrecht: Kluwer Academic Publisher, 1989.
2. Вильсон А.Дж. Энтропийные методы моделирования сложных систем. М.: Наука, 1978.
3. Николис Г., Пригожин И. Самоорганизация в неравновесных системах. М.: Мир, 1979.
4. Хакен Г. Информация и самоорганизация. Макроскопический подход к сложным системам. М.: УРСС, 2005.
5. Шредингер Э. Статистическая термодинамика. М.: ИЛ, 1948.
6. Крылов Н.С. Работы по обоснованию статистической физики. М. — Л.: Издательство АН СССР, 1950.
7. Хинчин А.Я. Математические основания статистической механики. М.-Ижевск: НИЦ «РХД», ИКИ, 2003.
8. Кац М. Вероятность и смежные вопросы в физике. М.: Мир, 1965.
9. Хуанг К. Статистическая механика. М.: Мир, 1966.
10. Рюэль Д. Статистическая механика. Строгие результаты. М.: Мир, 1971.
11. Корнфельд И.П., Синай Я.Г., Фомин С.В. Эргодическая теория. М.: Наука, 1980.
12. Evans L.C. Entropy and partial differential equations. Department of mathematics, UC Berkeley, 2003. http://math.berkeley.edu/~evans/
13. Минлос Р.А. Введение в математическую статистическую физику. М.: МЦНМО, 2002.
14. Козлов В.В. Ансамбли Гиббса и неравновесная статистическая механика. М.-Ижевск: НИЦ «РХД», ИКИ, 2008.
15. Марри Дж. Нелинейные дифференциальные уравнения в биологии. М.: Мир, 1983.
16. Свирежев Ю.М. Нелинейные волны, диссипативные структуры и катастрофы в экологии. М.: Наука, 1987.
17. Гардинер К.В. Стохастические методы в естественных науках. М.: Мир, 1986.
18. Веденяпин В.В. Кинетическая уравнения Больцмана и Власова. М.: Физматлит, 2001.
19. Малышев В.А., Пирогов С.А. Обратимость и необратимость в стохастической химической кинетике // УМН. 2008. Т. 63. № 1. С. 3-36.
20. Вайдлих В. Социодинамика: системный подход к математическому моделированию в социальных науках. М.: УРСС, 2010.
21. Castellano C, Fortunate S., Loreto V. Statistical physics of social behavior // Review of modern physics. 2009. V. 81. P. 591-646. arXiv:0710.3256v2
22. Занг В.-Б. Синергетическая экономика: время и перемены в нелинейной экономической теории. М.: Мир, 1999.
23. Dragulescu A., Yakovenko V.M. Statistical mechanics of money // The European Physical Journal B. 2000. V. 17. P. 723-729. arXiv: cond-mat /0001432 v4
24. Baldi P., Frasconi P., Smyth P. Modeling the Internet and the Web: Probabilistic methods and algorithms. Published by John Wiley & Sons, 2003.
25. Розоноэр Л.И. Обмен и распределение ресурсов (обобщенный термодинамический подход) I, II, III // Автоматика и телемеханика. 1973. № 5, № 6, № 8.
26. Горбань А.Н. Обход равновесия. Новосибирск: Наука, 1984.
27. Опойцев В.И. Нелинейная системостатика. М.: Наука, 1986.
28. Малишевский А.В. Качественные модели в теории сложных систем. М.: Наука, 1998.
29. Сергеев В.М. Пределы рациональности. М.: Фазис, 1999.
30. Попков Ю.С. Теория макросистем: равновесные модели. М.: УРСС, 1999.
31. Цирлин А.М. Методы оптимизации в необратимой термодинамике и микроэкономике. М.: Физматлит, 2003.
32. Швецов В.И., Алиев А.С. Математическое моделирование загрузки транспортных сетей. М.: УРСС, 2003.
33. Маслов В.П. Квантовая экономика. М.: Наука, 2006.
34. Олемской А.И. Синергетика сложных систем: Феноменология и статистическая теория. М.: КРАСАНД, 2009.
35. Веретенников А.Ю. Параметрическое и непараметрическое оценивание для цепей Маркова. М.: Изд-во ЦПИ при механико-математическом факультете МГУ, 2000.
36. Боровков А.А. Эргодичность и устойчивость случайных процессов. М.: УРСС, 1999.
37. Булинский А.В., Ширяев А.Н. Теория случайных процессов. М.: Физматлит; Лаборатория базовых знаний, 2003.
38. Кельберт М.Я., Сухов Ю.М. Вероятность и статистика в примерах и задачах. Т. 2: Марковские процессы как отправная точка теории случайных процессов и их приложения. М.: МЦНМО, 2010.
39. Вишневский В.М. Теоретические основы проектирования компьютерных сетей. М.: Техносфера, 2003.
40. Ивницкий В.А. Теория сетей массового обслуживания. М.: Физматлит, 2004.
41. The maximum entropy formalism, ed. by R.D. Levin, M. Tribus. Conf. Mass. Inst. Tech., Cambridge 1978. MIT Press, 1979.
42. International workshops on Bayesian inference and maximum entropy methods in science and engineering. AIP Conf. Proceedings (holds every year from 1980).
43. Kapur J.N. Maximum — entropy models in science and engineering. John Wiley & Sons, Inc., 1989.
44. Golan A., Judge G, Miller D. Maximum
entropy econometrics: Robust estimation with
limited data. Chichester, Wiley, 1996.
45. Fang S.-C, Rajasekera J.R., Tsao H.-S.J. Entropy optimization and mathematical programming. Kluwer’s International Series, 1997.
46. Богданов К.Ю. Прогулки с физикой. Библиотечка «Квант» В. 98. М.: Бюро Квантум, 2006 (глава 18).
47. Зорич В.А. Математический анализ задач естествознания. М.: МЦНМО, 2008.
48. Diaconis P. The Markov chain Monte Carlo revolution // Bulletin (New Series) of the AMS. 2009. V. 49, № 2. — P. 179-205. http://www. ams.org/journals/bull/2009-46-02/S0273-0979-08-01 238-X/S0273-0979-08-01238-X.pdf
49. Joulin A., Ollivier Y. Curvature, concentration and error estimates for Markov chain Monte Carlo // Ann. Prob. 2010. V. 38, № 6. — P. 2418-2442.
50. Кингман Дж. Пуассоновские процессы / под ред. А.М. Вершика. М.: МЦНМО, 2007.
51. Магарил-Ильяев Г.Г., Тихомиров В.М. Выпуклый анализ и его приложения. М.: УРСС, 2003.
52. Гасникова Е.В. Двойственные мультипликативные алгоритмы для задач энтропийно-линейного программирования // ЖВМ и МФ. 2009. Т. 49, № 3. — С. 453-464.
53. Нестеров Ю.Е. Введение в выпуклую оптимизацию. М.: МЦНМО, 2010.
54. Введение в математическое моделирование транспортных потоков: учеб. пособие / Гас-ников А.В., Кленов С.Л., Нурминский Е.А., Холодов Я.А., Шамрай Н.Б; Приложения: Бланк М.Л., Гасникова Е.В., Замятин А.А. и Малышев В.А., Колесников А.В., Райгородский А.М; под ред. А.В. Гасникова. — М.: МФТИ, 2010.
55. Flajolet P., Sedgewick R. Analytic combinatorics. Cambridge University Press, 2008. http://algo.inria.fr/flajolet/Publications/book.pdf
56. Вершик А.М., Шмидт А.А. Предельные меры, возникающие в асимптотической теории симметрических групп // ТВП. 1977. Т. 22. № 1. С. 72-88; 1978. Т. 23. № 1. С. 42-54.
57. Синай Я.Г. Вероятностный подход к анализу статистики выпуклых ломаных // Функц. анализ и его прил. 1994. Т. 28. № 2. С. 41-48.
5S. Колчин В.Ф. Случайные графы. М.: Физ-матлит, 2004.
59. Ledoux M. Concentration of measure phenomenon. Providence, RI, Amer. Math. Soc., 2001 (Math. Surveys Monogr. V. 89).
60. Алон Н., Спенсер Дж. Вероятностный метод. М.: Бином, 2007.
61. Якымив А.Л. Вероятностные приложения тауберовых теорем. М.: Наука, 2005.
62. Ширяев А.Н. Вероятность 1, 2. М.: МЦН-МО, 2007.
63. Malyshev V.A., Pirogov S.A., Rubco A.N. Random walks and chemical networks // Mosc. Math. J. 2004. V. 4, № 2. — P. 441-453.
64. Батищева Я.Г., Веденяпин В.В. II-й закон термодинамики для химической кинетики // Матем. моделирование. 2005. Т. 17, № 8. — С. 106-110.
65. Веденяпин В.В., Орлов Ю.Н. О законах сохранения для полиномиальных гамильтонианов и для дискретных моделей уравнения Больцмана // ТМФ. 1999. Т. 121, № 2. — С. 307-315.
66. Рыбко А.Н. Пуассоновская гипотеза для больших симметричных коммуникационных сетей // Глобус. Общематематический семинар / под ред. М.А. Цфасмана и В.В. Прасолова. № 4. М.: МЦНМО, 2009. — С. 105-126.
67. Федорюк М.В. Метод перевала. М.: УРСС, 2010.
Поступила в редакцию 15.10.2010.
