Журнал вычислительной математики и математической физики, 2019, T. 59, № 1, стр. 21-36

Универсальный метод поиска равновесий и стохастических равновесий в транспортных сетях
Д. Р. Баймурзина, А. В. Гасников, Е. В. Гасникова, П. Е. Двуреченский, Е. И. Ершов, М. Б. Кубентаева, А. А. Лагуновская

Д. Р. Баймурзина^1, 2, *, А. В. Гасников^1, 3, **, Е. В. Гасникова^1, ***, П. Е. Двуреченский^3, 4, ****, Е. И. Ершов^3, *****, М. Б. Кубентаева^1, ******, А. А. Лагуновская^1, *******

¹ МФТИ
141700 М.о., Долгопрудный, Институтский пер., 9, Россия

² “Сколково”
143026 Москва, ул. Нобеля, 3, Россия

³ ИППИ РАН
127051 Москва, Б. Каретный пер., 19, стр. 1, Россия

⁴ WIAAS
10117 Berlin, Mohrenstr, 39, Germany

^* E-mail: dilyara.rimovna@gmail.com
^** E-mail: gasnikov.av@mipt.ru
^*** E-mail: egasnikova@yandex.ru
^**** E-mail: dvurechensky@iitp.ru
^***** E-mail: e.i.ershov@gmail.com
^****** E-mail: kubikmeruza@yandex.ru
^******* E-mail: a.lagunovskaya@phystech.edu

Поступила в редакцию 19.01.2017
После доработки 04.12.2017

DOI: 10.1134/S0044466919010022

Полный текст (PDF)

Аннотация

Предложен универсальный способ поиска обычных и стохастических равновесий в популяционных играх загрузки. Рассматриваются модели равновесного распределения потоков по путям Бэкмана и стабильной динамики. Поиск стохастических равновесий Нэша(–Вардропа) приводит к решению энтропийно регуляризованных задач выпуклой оптимизации. Данная работа посвящена поиску эффективного решения такого рода задачи, а точнее двойственных к ним, с помощью недавно предложенного прямо-двойственного универсального градиентного метода, оптимально и адаптивно настраивающегося на гладкость решаемой задачи. Библ. 45. Фиг. 2. Табл. 1.

Ключевые слова: транспортные потоки, универсальный метод подобных треугольников, двойственная задача, модель Бэкмана, модель стабильной динамики.

1. ВВЕДЕНИЕ

В [1] было анонсировано, что в цикле последующих статей планируется описать новые эффективные численные методы поиска равновесий (и стохастических равновесий [2]) в популярных моделях распределения транспортных потоков по путям (ребрам) графа транспортной сети [2]–[5]. В частности, это планировалось сделать для модели Бэкмана [3], [5] (называемой также BMW-моделью) и модели стабильной динамики [4], [5] (называемой также моделью Нестерова–де Пальмы). В цикле работ [6]–[14] этот план начал осуществляться. Однако на данный момент остался открытым главный вопрос: какие из предложенных методов являются наилучшими? Одним из основных результатов данной работы является попытка сопоставления и унификации имеющихся здесь результатов с целью получения ответа на поставленный вопрос. Оказалось, что для поиска обычного (не стохастического) равновесия в модели Бэкмана наилучшим методом (по имеющимся на данный момент представлениям) является классический метод условного градиента Франк–Вульфа (см. [3], [6]), а для всех остальных случаев (трех: поиск стохастических равновесий в модели Бэкмана и поиск обычных и стохастических равновесий в модели стабильной динамики) разумно строить двойственную задачу (эта задача оказывается негладкой при поиске обычных (не стохастических) равновесий) и решать ее универсальным градиентным методом Ю.Е. Нестерова (см. [15] и также более поздние версии [16], [17]). Затем, пользуясь прямодвойственностью этого метода (см. [17]), можно восстанавливать решение исходной (прямой задачи). Интересно отметить, что численные эксперименты (см. п. 7), проведенные с универсальным градиентным методом для поиска обычных (не стохастических) равновесий показали, что время работы метода пропорционально $ \sim {\kern 1pt} {{\tilde {\varepsilon }}^{{ - 1}}}$, где $\tilde {\varepsilon }$ – желаемая относительная точность (по функции) решения задачи, а не $ \sim {\kern 1pt} {{\tilde {\varepsilon }}^{{ - 2}}}$, как это можно было ожидать, исходя из теории решения негладких задач выпуклой оптимизации [18]. В отличие от специальных примеров, собранных в [15], [16], на которых наблюдался аналогичный эффект, в данной статье обнаружен, по-видимому, первый совершенно реальный такой пример с существенно негладкой функцией, имеющей огромное число всевозможных изломов.

В разд. 2 приводятся экстремальные принципы, описывающие равновесную конфигурацию в транспортных сетях в четырех рассматриваемых случаях. Эти принципы формулируются как сепарабельные задачи выпуклой оптимизации (с композитом вида энтропии в случае поиска стохастических равновесий) при аффинных ограничениях. В основном в изложении используются результаты работ [1], [5], [8], [9], [11], [19]–[26].

В разд. 3 с помощью достаточно стандартной техники выпуклого анализа строятся двойственные задачи, чтобы не проектироваться на аффинные ограничения при решении прямых задач (последняя операция является дорогой [18]). Во многом здесь используются результаты работ [8], [9], [12].

В разд. 4, следуя [17], излагается специальный вариант универсального градиентного метода – универсальный метод подобных треугольников (УМПТ), в котором (в отличие от [15], [16]) используется только одна проекция на каждой итерации. Особое внимание уделяется прямодвойственности этого метода. Это нашло отражение в необычной форме записи того, как сходится метод.

В разд. 5 УМПТ применяется для решения двойственных задач, описанных в разд. 3. Приводятся формулы восстановления решения прямой задачи. В основном изложение следует [9], [12]. Формулы восстановления решения задачи поиска равновесия в модели стабильной динамики, исходя из последовательности, сгенерированной УМПТ при решении двойственной задачи, по-видимому, приводятся впервые.

В разд. 5 внимание было сконцентрировано на том, как восстанавливать решения прямых задач, т.е. задач из разд. 2, исходя из приближенного решения двойственных задач, т.е. задач из разд. 3. В разд. 6 результаты разд. 5 пополняются исследованием того, сколько итераций будет делать УМПТ (оценка констант метода) для достижения желаемой точности решения (прямой и двойственной задачи одновременно), и какая при этом будет стоимость одной итерации. Последняя задача сводится (при поиске стохастических равновесий) к вычислению значений и градиентов характеристических функций на транспортном графе [5], [9], [27] (для вычисления градиентов используется теория быстрого автоматического дифференцирования [28]) или к поиску кратчайших путей в транспортном графе [29] (при поиске обычных равновесий).

В разд. 7 результаты всех предыдущих разделов собираются вместе. Формулируется теорема, в которой описаны полученные в данной работе (верхние) теоретические оценки скорости сходимости УМПТ во всех четырех случаях. Приводятся результаты численных экспериментов. Эти результаты сопоставляются с ранее известными.

2. ЭКСТРЕМАЛЬНЫЕ ПРИНЦИПЫ ДЛЯ ПОИСКА РАВНОВЕСИЙ В ТРАНСПОРТНЫХ СЕТЯХ

Рассмотрим транспортную сеть, которую будем представлять ориентированным графом $\left\langle {V,E} \right\rangle $, где $V$ – множество вершин (как правило, можно считать, что ${{\left| E \right|} \mathord{\left/ {\vphantom {{\left| E \right|} 4}} \right. \kern-0em} 4} \leqslant \left| V \right| \leqslant \left| E \right|$), а $E$ – множество ребер, $\left| E \right| = n$. Обозначим множество пар источник-сток $w = \left( {i,j} \right)$ через $OD$, через ${{d}_{w}}$ – корреспонденцию, отвечающую паре $w$, ${{x}_{p}}$ – поток по пути $p$; ${{P}_{w}}$ – множество путей, отвечающих корреспонденции $w$ (начинающихся в $i$ и заканчивающихся в $j$), $P = \bigcup\nolimits_{w \in OD}^{} {{{P}_{w}}} $ – множество всех путей. Затраты на прохождения ребра $e \in E$ описываются функцией ${{\tau }_{e}}\left( {{{f}_{e}}} \right)$, где ${{f}_{e}}$ – поток по ребру $e$.

Опишем марковскую логит-динамику (также говорят гиббсовскую динамику) в повторяющейся игре загрузки графа транспортной сети [19]. Пусть каждой корреспонденции отвечает ${{d}_{w}}M$ агентов ($M \gg 1$),${{\tau }_{e}}\left( {{{f}_{e}}} \right): = {{\tau }_{e}}\left( {{{{{f}_{e}}} \mathord{\left/ {\vphantom {{{{f}_{e}}} M}} \right. \kern-0em} M}} \right)$. Пусть имеется $TN$ шагов ($N \gg 1$). Каждый агент независимо от остальных на шаге $t + 1$ выбирает с вероятностью

$\frac{\lambda }{N}\frac{{\exp \left( { - {{G_{p}^{t}} \mathord{\left/ {\vphantom {{G_{p}^{t}} \gamma }} \right. \kern-0em} \gamma }} \right)}}{{\sum\limits_{q \in {{P}_{w}}} {\exp \left( { - {{G_{q}^{t}} \mathord{\left/ {\vphantom {{G_{q}^{t}} \gamma }} \right. \kern-0em} \gamma }} \right)} }}.$

Путь $p \in {{P}_{w}}$, где $G_{p}^{t}$ – затраты на пути $p$ на шаге $t$ ($G_{p}^{0} \equiv 0$), а с вероятностю $1 - {\lambda \mathord{\left/ {\vphantom {\lambda N}} \right. \kern-0em} N}$ путь, который использовался на шаге $t$. Такая динамика отражает ограниченную рациональность агентов (см. разд. 3), и часто используется в популяционной теории игр [19] и теории дискретного выбора [20]. Оказывается, эта марковская динамика в пределе $N \to \infty $ превращается в марковскую динамику в непрерывном времени (вырождающуюся при $\gamma \to 0{\kern 1pt} + $ в динамику наилучших ответов [19]), которая в свою очередь допускает два предельных перехода (обоснование перестановочности этих пределов см. в [21]): $T \to \infty $, $M \to \infty $ или $M \to \infty $, $T \to \infty $. При первом порядке переходов мы сначала ($T \to \infty $) согласно эргодической теореме для марковских процессов (в нашем случае марковский процесс – модель стохастической химической кинетики с унарными реакциями в условиях детального баланса [22]–[24]) приходим к финальной (=стационарной) вероятностной мере, имеющей в основе мультиномиальное распределение. С ростом числа агентов ($M \to \infty $) эта мера концентрируется около наиболее вероятного состояния, поиск которого сводится к решению задачи (1) ниже. Функционал в этой задаче оптимизации с точностью до потенцирования и мультипликативных и аддитивных констант соответствует исследуемой стационарной мере – т.е. это функционал Санова [24], [25]. При обратном порядке переходов мы сначала осуществляем так называемый канонический скейлинг [21], [22], приводящий к детерминированной кинетической динамике, описываемой СОДУ на $x$, а затем ($T \to \infty $) ищем аттрактор получившейся СОДУ. Глобальным аттрактором оказывается неподвижная точка, которая определяется решением задачи (1) (cм. ниже). Более того, функционал, стоящий в (1), является функцией Ляпунова полученной кинетической динамики (т.е. функционалом Больцмана). Последнее утверждение – достаточно общий факт (функционал Санова, является функционалом Больцмана), верный при намного более общих условиях(см. [24]).

Итак, рассматривается следующая задача поиска стохастического равновесия Нэша–Вардропа в модели Бэкмана равновесного распределения транспортных потоков по путям (см., например, [1], [2]):

(1)

$\sum\limits_{e \in E} {{{\sigma }_{e}}\left( {{{f}_{e}}} \right)} + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {{{x}_{p}}\ln \left( {{{{{x}_{p}}} \mathord{\left/ {\vphantom {{{{x}_{p}}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } \to \mathop {\min }\limits_{f = \Theta x,\;x \in X} ,$

где $\gamma > 0$; ${{\sigma }_{e}}\left( {{{f}_{e}}} \right) = \int_0^{{{f}_{e}}} {{{\tau }_{e}}\left( z \right)dz} $ – выпуклые функции;

$\Theta = {{\left\| {{{\delta }_{{ep}}}} \right\|}_{{e \in E,p \in P}}} = {{\left\| {{{\Theta }^{{\left\langle p \right\rangle }}}} \right\|}_{{p \in P}}},\quad {{\delta }_{{ep}}} = \left\{ \begin{gathered} 1,\quad e \in p, \hfill \\ 0,\quad e \notin p, \hfill \\ \end{gathered} \right.$

$X = \left\{ {x \geqslant 0\,:\;\;\sum\limits_{p \in {{P}_{w}}} {{{x}_{p}}} = {{d}_{w}},\;w \in OD} \right\}$

есть прямое произведение симплексов;

в качестве ${{\tau }_{e}}\left( {{{f}_{e}}} \right)$ обычно выбирают BPR-функции (см. [1], [3], [5], [30]):

${{\tau }_{e}}\left( {{{f}_{e}}} \right) = {{\bar {t}}_{e}}\left( {1 + \rho {{{\left( {{{{{f}_{e}}} \mathord{\left/ {\vphantom {{{{f}_{e}}} {{{{\bar {f}}}_{e}}}}} \right. \kern-0em} {{{{\bar {f}}}_{e}}}}} \right)}}^{4}}} \right).$

В пределе модели стабильной динамики [1] имеем

$\tau _{e}^{\mu }\left( {{{f}_{e}}} \right) = {{\bar {t}}_{e}}\left( {1 + \rho {{{\left( {{{{{f}_{e}}} \mathord{\left/ {\vphantom {{{{f}_{e}}} {{{{\bar {f}}}_{e}}}}} \right. \kern-0em} {{{{\bar {f}}}_{e}}}}} \right)}}^{{{1 \mathord{\left/ {\vphantom {1 \mu }} \right. \kern-0em} \mu }}}}} \right)\;\xrightarrow[{\mu \to 0 + }]{}\;\left\{ \begin{gathered} {{{\bar {t}}}_{e}},\quad 0 \leqslant {{f}_{e}} < {{{\bar {f}}}_{e}}, \hfill \\ \left[ {{{{\bar {t}}}_{e}},\infty } \right),\quad {{f}_{e}} = {{{\bar {f}}}_{e}}, \hfill \\ \end{gathered} \right.$

${{d\tau _{e}^{\mu }\left( {{{f}_{e}}} \right)} \mathord{\left/ {\vphantom {{d\tau _{e}^{\mu }\left( {{{f}_{e}}} \right)} {d{{f}_{e}}}}} \right. \kern-0em} {d{{f}_{e}}}}\;\xrightarrow[{\mu \to 0 + }]{}\;0,\quad 0 \leqslant {{f}_{e}} < {{\bar {f}}_{e}},$

задача перепишется в виде (см. [1])

(2)

3. ПЕРЕХОД К ДВОЙСТВЕННОЙ ЗАДАЧЕ

Запишем двойственную задачу к (1) [1], [6], [9] (далее мы используем обозначение ${\text{dom}}\sigma {\text{*}}$ – область определения сопряженной к $\sigma $ функции):

(3)

где (oбратим внимание, что вектор распределения потоков по путям $x$ при поиске стохастического равновесия ($\gamma > 0$) получается не разреженным в отличие от поиска обычного равновесия Нэша(–Вардропа) [6] ($\gamma \to 0 + $). Как следствие, чтобы вычислить этот вектор, требуются затраты, существенно зависящие от потенциально огромной размерности $\left| P \right|$. К счастью, в приложениях, как правило, не требуется знание этого вектора, достаточно определить вектор потоков на ребрах $f$, который, как мы увидим ниже, может быть вычислен намного эффективнее (в частности, с затратами, не зависящими от $\left| P \right|$)

(4)

$\begin{gathered} \psi \left( t \right) = \sum\limits_{w \in OD} {{{d}_{w}}{{\psi }_{w}}\left( t \right)} ,\quad {{\psi }_{w}}\left( t \right) = \ln \left( {\sum\limits_{p \in {{P}_{w}}} {\exp \left( { - \sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \right)} } \right), \\ f = - \nabla \gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right),\quad {{x}_{p}} = {{d}_{w}}\frac{{\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \right)}}{{\sum\limits_{q \in {{P}_{w}}} {\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{eq}}}{{t}_{e}}} } \right)} }},\quad p \in {{P}_{w}}, \\ \end{gathered} $

для

${{\tau }_{e}}\left( {{{f}_{e}}} \right) = {{\bar {t}}_{e}}\left( {1 + \rho {{{\left( {\frac{{{{f}_{e}}}}{{{{{\bar {f}}}_{e}}}}} \right)}}^{{\frac{1}{\mu }}}}} \right),$

имеем [1] (BPR-функция, см. разд. 2, получается при $\mu = {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0em} 4}$)

$\sigma _{e}^{*}\left( {{{t}_{e}}} \right) = \mathop {\sup }\limits_{{{f}_{e}} \geqslant 0} \left( {\left( {{{t}_{e}} - {{{\bar {t}}}_{e}}} \right){{f}_{e}} - {{{\bar {t}}}_{e}}\frac{\mu }{{1 + \mu }}\rho \frac{{{{f}_{e}}^{{1 + \frac{1}{\mu }}}}}{{{{{\bar {f}}}_{e}}^{{\frac{1}{\mu }}}}}} \right) = {{\bar {f}}_{e}}{{\left( {\frac{{{{t}_{e}} - {{{\bar {t}}}_{e}}}}{{{{{\bar {t}}}_{e}}\rho }}} \right)}^{\mu }}\frac{{\left( {{{t}_{e}} - {{{\bar {t}}}_{e}}} \right)}}{{1 + \mu }}.$

Собственно, формула (3) есть не что иное, как отражение формулы $f = - \nabla \gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ и связи ${{t}_{e}} = {{\tau }_{e}}\left( {{{f}_{e}}} \right)$, $e \in E$. Действительно, по формуле Демьянова–Данскина–Рубинова (см. [31]) имеем

$\frac{{d\sigma _{e}^{*}\left( {{{t}_{e}}} \right)}}{{d{{t}_{e}}}} = \frac{d}{{d{{t}_{e}}}}\mathop {\max }\limits_{{{f}_{e}} \geqslant 0} \left\{ {{{t}_{e}}{{f}_{e}} - \int\limits_0^{{{f}_{e}}} {{{\tau }_{e}}\left( z \right)dz} } \right\} = {{f}_{e}}\,{\text{:}}\;\;{{t}_{e}} = {{\tau }_{e}}\left( {{{f}_{e}}} \right).$

В свою очередь, формула $f = - \nabla \gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ может интерпретироваться как следствие соотношений $f = \Theta x$ и формулы распределения Гиббса (логит-распределения)

${{x}_{p}} = {{d}_{w}}\frac{{\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \right)}}{{\sum\limits_{q \in {{P}_{w}}} {\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{eq}}}{{t}_{e}}} } \right)} }},\quad p \in {{P}_{w}}.$

При такой интерпретации связь задачи (3) с логит-динамикой, порождающей стохастические равновесия, наиболее наглядна. Последняя формула, в виду того, что ${{g}_{p}}\left( t \right) = \sum\nolimits_{e \in E}^{} {{{\delta }_{{ep}}}{{t}_{e}}} $ – затраты на пути $p$ на графе $\left\langle {V,E} \right\rangle $, ребра которого взвешены, $t$ – есть отражение следующего принципа поведения (ограниченной рациональности агентов [20]): каждый агент $k$ (пользователь транспортной сети), отвечающий корреспонденции $w \in OD$, выбирает маршрут следования $p \in {{P}_{w}}$, если

$p = \arg \mathop {\max }\limits_{q \in {{P}_{w}}} \left\{ { - {{g}_{q}}\left( t \right) + \xi _{q}^{k}} \right\},$

где независимые случайные величины $\xi _{q}^{k}$, имеют одинаковое двойное экспоненциальное распределение, также называемое распределением Гумбеля (см. [19], [20]):

$P\left( {\xi _{q}^{k} < \zeta } \right) = \exp \left\{ { - {{e}^{{ - {\zeta \mathord{\left/ {\vphantom {\zeta \gamma }} \right. \kern-0em} \gamma } - E}}}} \right\},\quad \gamma > 0,$

распределение Гумбеля можно объяснить исходя из идемпотентного аналога центральной предельной теоремы (вместо суммы случайных величин – максимум) для независимых случайных величин с экспоненциальным и более быстро убывающим правым хвостом [32]. Распределение Гумбеля возникает в данном контексте, например, если при принятии решения водитель собирает информацию с большого числа разных (независимых) зашумленных источников, ориентируясь на худшие прогнозы по каждому из путей).

Отметим также, что если взять $E \approx 0.5772$ – константа Эйлера, то ${\text{M}}\left[ {\xi _{q}^{k}} \right] = 0$, $D\left[ {\xi _{q}^{k}} \right] = {{{{\gamma }^{2}}{{\pi }^{2}}} \mathord{\left/ {\vphantom {{{{\gamma }^{2}}{{\pi }^{2}}} 6}} \right. \kern-0em} 6}$. Распределение Гиббса получается в пределе, когда число агентов на каждой корреспонденции стремится к бесконечности (случайность исчезает и описание переходит на средние величины).

Сформулируем главный вывод из написанного выше. Если каждый пользователь сориентирован на вектор затрат $t$ на ребрах $E$, одинаковый для всех пользователей, и пытается выбрать кратчайший путь, исходя из зашумленной информации, то такое поведение пользователей в пределе, когда их число стремится к бесконечности $M \to \infty $, приводит к описанию распределения пользователей по путям/ребрам (4). Равновесная конфигурация, описываемая решением задачи (3), характеризуется тем, что по вектору $t$ вычисляется согласно формуле (4) такой вектор $f = \Theta x$, что имеет место соотношение $t = {{\left\{ {{{\tau }_{e}}\left( {{{f}_{e}}} \right)} \right\}}_{{e \in E}}}$.

Заметим, что при $\gamma \to 0{\kern 1pt} + $ распределение водителей по путям вырождается, и все водители (агенты) будут использовать только кратчайшие пути

(5)

$ - \mathop {\lim }\limits_{\gamma \to 0 + } \gamma {{\psi }_{w}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = \mathop {\min }\limits_{p \in {{P}_{w}}} {{g}_{p}}\left( t \right),$

(6)

$ - \mathop {\lim }\limits_{\gamma \to 0 + } \nabla \gamma {{\psi }_{w}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = {\text{conv}}\left\{ {{{{\left( {0,1,\;...,\;0,1} \right)}}^{{\text{T}}}}\;,...,\;{{{(1,1,\;...,\;0,0)}}^{{\text{T}}}}} \right\},$

в правой части формулы (6) стоит субградиент (негладкой) функции длины кратчайшего пути, отвечающего корреспонденции $w$. Этот субградиент (в общем случае) есть множество, представимое в виде выпуклой оболочки векторов, отвечающих всевозможным кратчайшим путям, начинающимся в вершине $i$ и заканчивающимся в вершине $j$, где $w = \left( {i,j} \right)$, в рассматриваемом транспортном графе, ребра которого взвешены вектором $t$. Единицы стоят в векторах на местах, которые отвечают ребрам, входящим в данный кратчайший путь. Если кратчайший путь единственен, то субградиент превращается в обычный градиент. Далее (в УМПТ) можно выбирать произвольный элемент субградиента. От этого выбора приводимые ниже оценки не будут зависеть.

Полезно также иметь в виду, что [19], [20]

$\gamma {{\psi }_{w}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = {{{\text{M}}}_{{{{{\left\{ {{{\xi }_{p}}} \right\}}}_{{p \in {{P}_{w}}}}}}}}\left[ {\mathop {\max }\limits_{p \in {{P}_{w}}} \left\{ { - {{g}_{p}}\left( t \right) + {{\xi }_{p}}} \right\}} \right].$

В пределе стабильной динамики $\mu \to 0{\kern 1pt} + $ задача (3) вырождается в задачу

$\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\bar {f},t - \bar {t}} \right\rangle \to \mathop {\min }\limits_{t \geqslant \bar {t}} .$

В пределе $\gamma \to 0{\kern 1pt} + $ все приведенные выше формулы переходят в соответствующие формулы для классических (не стохастических) моделей Бэкмана и стабильной динамики, см., например, [6].

4. УНИВЕРСАЛЬНЫЙ МЕТОД ПОДОБНЫХ ТРЕУГОЛЬНИКОВ

Ниже, следуя [17], изложен возможный способ решения двойственной задачи (3), позволяющий при этом восстанавливать решение прямой задачи ((1) или (2)), исходя из подхода, описанного в [33], [34] (см. также первоисточники [35], [36]). Рассмотрим общую задачу выпуклой композитной оптимизации на множестве простой структуры

(7)

$F\left( t \right) = \underbrace {\Phi \left( t \right)}_{\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)} + \underbrace {h\left( t \right)}_{\sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \to \mathop {\min }\limits_{t \in Q} .$

Положим ${{R}^{2}} = \frac{1}{2}\left\| {{{t}_{*}} - {{y}^{0}}} \right\|_{2}^{2}$, где ${{y}^{0}} = \bar {t}$, ${{t}_{*}}$ – решение задачи (7) (если решение не единственно, то выбирается то, которое доставляет минимум $\left\| {{{t}_{*}} - {{y}^{0}}} \right\|_{2}^{2}$).

Пусть (${{L}_{\nu }} \leqslant \infty $, т.е. допускается равенство бесконечности)

(8)

${{\left\| {\nabla \Phi \left( t \right) - \nabla \Phi \left( y \right)} \right\|}_{2}} \leqslant {{L}_{\nu }}\left\| {t - y} \right\|_{2}^{\nu },\quad \nu \in \left[ {0,1} \right],\quad {{L}_{0}} < \infty .$

Положим

$\begin{gathered} {{\varphi }_{0}}\left( t \right) = {{\alpha }_{0}}\left[ {\Phi \left( {{{y}^{0}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{0}}} \right),t - {{y}^{0}}} \right\rangle + h\left( t \right)} \right] + \frac{1}{2}\left\| {t - {{y}^{0}}} \right\|_{2}^{2}, \\ {{\varphi }_{{k + 1}}}\left( t \right) = {{\varphi }_{k}}\left( t \right) + {{\alpha }_{{k + 1}}}\left[ {\Phi \left( {{{y}^{{k + 1}}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{{k + 1}}}} \right),t - {{y}^{{k + 1}}}} \right\rangle + h\left( t \right)} \right]. \\ \end{gathered} $

Начнем описание УМПТ с самой первой итерации. Положим

${{A}_{0}} = {{\alpha }_{0}} = {1 \mathord{\left/ {\vphantom {1 {L_{0}^{0}}}} \right. \kern-0em} {L_{0}^{0}}},\quad k = 0,\quad {{j}_{0}} = 0;\quad {{t}^{0}} = {{u}^{0}} = \arg \mathop {\min }\limits_{t \in Q} {{\varphi }_{0}}\left( t \right).$

До тех пор пока имеем

$\Phi \left( {{{t}^{0}}} \right) > \Phi \left( {{{y}^{0}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{0}}} \right),{{t}^{0}} - {{y}^{0}}} \right\rangle + \frac{{L_{0}^{{{{j}_{0}}}}}}{2}\left\| {{{t}^{0}} - {{y}^{0}}} \right\|_{2}^{2} + \frac{{{{\alpha }_{0}}}}{{2{{A}_{0}}}}\varepsilon ,$

где

${{t}^{0}}: = {{u}^{0}}: = \arg \mathop {\min }\limits_{t \in Q} {{\varphi }_{0}}\left( t \right),\quad \left( {{{A}_{0}}: = } \right){{\alpha }_{0}}: = \frac{1}{{L_{0}^{{{{j}_{0}}}}}},$

выполнять

${{j}_{0}}: = {{j}_{0}} + 1;\quad L_{0}^{{{{j}_{0}}}}: = {{2}^{{{{j}_{0}}}}}L_{0}^{0}.$

Универсальный метод подобных треугольников

Шаг 1. $L_{{k + 1}}^{0} = {{L_{k}^{{{{j}_{k}}}}} \mathord{\left/ {\vphantom {{L_{k}^{{{{j}_{k}}}}} 2}} \right. \kern-0em} 2}$, ${{j}_{{k + 1}}} = 0$.

Шаг 2. $\left\{ \begin{gathered} {{\alpha }_{{k + 1}}}: = \frac{1}{{2L_{{k + 1}}^{{{{j}_{{k + 1}}}}}}} + \sqrt {\frac{1}{{4{{{\left( {L_{{k + 1}}^{{{{j}_{{k + 1}}}}}} \right)}}^{2}}}} + \frac{{{{A}_{k}}}}{{L_{{k + 1}}^{{{{j}_{{k + 1}}}}}}}} ,\quad {{A}_{{k + 1}}}: = {{A}_{k}} + {{\alpha }_{{k + 1}}}; \hfill \\ {{y}^{{k + 1}}}: = \frac{{{{\alpha }_{{k + 1}}}{{u}^{k}} + {{A}_{k}}{{t}^{k}}}}{{{{A}_{{k + 1}}}}},\quad {{u}^{{k + 1}}}: = \arg \mathop {\min }\limits_{t \in Q} {{\varphi }_{{k + 1}}}\left( t \right),\quad {{t}^{{k + 1}}}: = \frac{{{{\alpha }_{{k + 1}}}{{u}^{{k + 1}}} + {{A}_{k}}{{t}^{k}}}}{{{{A}_{{k + 1}}}}}. \hfill \\ \end{gathered} \right.$

До тех пор пока

$\Phi \left( {{{y}^{{k + 1}}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{{k + 1}}}} \right),{{t}^{{k + 1}}} - {{y}^{{k + 1}}}} \right\rangle + \frac{{L_{{k + 1}}^{{{{j}_{{k + 1}}}}}}}{2}\left\| {{{t}^{{k + 1}}} - {{y}^{{k + 1}}}} \right\|_{2}^{2} + \frac{{{{\alpha }_{{k + 1}}}}}{{2{{A}_{{k + 1}}}}}\varepsilon < \Phi \left( {{{t}^{{k + 1}}}} \right),$

выполнять

${{j}_{{k + 1}}}: = {{j}_{{k + 1}}} + 1;\quad L_{{k + 1}}^{{{{j}_{{k + 1}}}}} = 2L_{{k + 1}}^{{{{j}_{k}}}};$

Шаг 3. $k: = k + 1$ и go to 1.

УМПТ сходится согласно оценке

(9)

${{A}_{N}}F\left( {{{t}^{N}}} \right) \leqslant \mathop {\min }\limits_{t \in Q} \left\{ {\frac{1}{2}\left\| {t - {{y}^{0}}} \right\|_{2}^{2} + \sum\limits_{k = 0}^N {{{\alpha }_{k}}\left[ {\Phi \left( {{{y}^{k}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{k}}} \right),t - {{y}^{k}}} \right\rangle + h\left( t \right)} \right]} } \right\}.$

Отсюда

$\begin{gathered} {{A}_{N}}F\left( {{{t}^{N}}} \right) \leqslant \mathop {\min }\limits_{t \in Q} \left\{ {\frac{1}{2}\left\| {t - {{y}^{0}}} \right\|_{2}^{2} + \sum\limits_{k = 0}^N {{{\alpha }_{k}}\left[ {\Phi \left( {{{y}^{k}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{k}}} \right),t - {{y}^{k}}} \right\rangle + h\left( t \right)} \right]} } \right\} + \frac{\varepsilon }{2} \leqslant \\ \leqslant \;\frac{1}{2}\left\| {{{t}_{*}} - {{y}^{0}}} \right\|_{2}^{2} + \sum\limits_{k = 0}^N {{{\alpha }_{k}}\underbrace {\left[ {\Phi \left( {{{y}^{k}}} \right) + \left\langle {\nabla \Phi \left( {{{y}^{k}}} \right),{{t}_{*}} - {{y}^{k}}} \right\rangle + h\left( {{{t}_{*}}} \right)} \right]}_{ \leqslant \Phi \left( {{{t}_{*}}} \right) + h\left( {{{t}_{*}}} \right)}} + \frac{\varepsilon }{2} \leqslant \\ \leqslant \;\frac{1}{2}\left\| {{{t}_{*}} - {{y}^{0}}} \right\|_{2}^{2} + \sum\limits_{k = 0}^N {{{\alpha }_{k}}F\left( {{{t}_{*}}} \right)} + \frac{\varepsilon }{2} = {{R}^{2}} + {{A}_{N}}F\left( {{{t}_{*}}} \right) + \frac{\varepsilon }{2}, \\ \end{gathered} $

т.е.

$F\left( {{{t}^{N}}} \right) - F\left( {{{t}_{*}}} \right) \leqslant \frac{{{{R}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}.$

В [17] доказано, что ${{A}_{N}} \simeq {{2{{R}^{2}}} \mathord{\left/ {\vphantom {{2{{R}^{2}}} \varepsilon }} \right. \kern-0em} \varepsilon }$ при

(10)

$N \simeq \mathop {\inf }\limits_{\nu \in \left[ {0,1} \right]} {{\left( {\frac{{{{L}_{\nu }}{{{\left( {16R} \right)}}^{{1 + \nu }}}}}{\varepsilon }} \right)}^{{\frac{2}{{1 + 3\nu }}}}},$

т.е. при таком $N$ справедливо неравенство

$F\left( {{{t}^{N}}} \right) - \mathop {\min }\limits_{t \in Q} F\left( t \right) \leqslant \varepsilon .$

При этом среднее число вычислений значения функции на одной итерации будет $ \approx {\kern 1pt} 4$, а градиента функции $ \approx {\kern 1pt} 2$.

5. ПРИЛОЖЕНИЕ УМПТ К ПОИСКУ РАВНОВЕСИЙ В ТРАНСПОРТНЫХ СЕТЯХ

Задачи, получающиеся из (3) при $\mu > 0$ и $\mu \to 0{\kern 1pt} + $, можно не различать по сложности при композитном подходе к ним (см. [17], [37]). В задаче, отвечающей $\mu \to 0{\kern 1pt} + $, сепарабельный композит $h\left( t \right)$ проще сепарабельного композита задачи с $\mu > 0$, но зато дополнительно добавляется сепарабельное ограничение простой структуры (для BPR-функций ($\mu = 0.25$) все сводится к решению уравнения четвертой степени, см. разд. 3 и [12]).

В любом случае основные затраты на каждой итерации при использовании УМПТ в композитном варианте [17], [37] связаны с необходимостью расчета градиента гладкой функции $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$. Тому, как это можно эффективно делать, будет посвящен разд. 6. Здесь же отметим, что константы Гёльдера (см. (8)) градиента этой функции в 2-норме могут быть оценены в виде (см. [38]–[40])

(11)

${{L}_{1}} \leqslant \frac{1}{\gamma }\sum\limits_{w \in OD} {{{d}_{w}}} \mathop {\max }\limits_{p \in {{P}_{w}}} \left\| {{{\Theta }^{{\left\langle p \right\rangle }}}} \right\|_{2}^{2} \leqslant \frac{{Hd}}{\gamma },$

(12)

${{L}_{0}} \leqslant 2\sum\limits_{w \in OD} {{{d}_{w}}} \mathop {\max }\limits_{p \in {{P}_{w}}} {{\left\| {{{\Theta }^{{\left\langle p \right\rangle }}}} \right\|}_{2}} \leqslant 2\sqrt H d\quad ({\text{п р и }}\;\gamma \to 0 + ),$

где $d = \sum\nolimits_{w \in OD}^{} {{{d}_{w}}} $, $H = \mathop {\max }\limits_{p \in P} \left\| {{{\Theta }^{{\left\langle p \right\rangle }}}} \right\|_{2}^{2}$, т.е. $H$ – максимальное число ребер в пути. Как правило, можно считать, что $H = {\rm O}\left( {\sqrt n } \right)$ – диаметр графа с манхетенской структурой, т.е. для квадратной решетки с $n$ ребрами.

В данном разделе внимание будет сосредоточено на том, как с помощью формулы (9) восстанавливать решение прямой задачи (см. (1) или (2)), исходя из последовательности, сгенерированной УМПТ при решении двойственной задачи (3). Описанные далее способы имеют много общего с методами, описанными в [12], [33], [34].

Положим

${{f}^{k}} = - \nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),\quad x_{p}^{k} = {{d}_{w}}\frac{{\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{ep}}}y_{e}^{k}} } \right)}}{{\sum\limits_{q \in {{P}_{w}}} {\exp \left( { - \frac{1}{\gamma }\sum\limits_{e \in E} {{{\delta }_{{eq}}}y_{e}^{k}} } \right)} }},\quad p \in {{P}_{w}},\quad w \in OD,$

${{\bar {f}}^{N}} = \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{a}_{k}}{{f}^{k}}} ,\quad {{\bar {x}}^{N}} = \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}{{x}^{k}}} .$

Пусть (минимум в (13) достигается во внутренней точке множества $\operatorname{dom} \sigma {\text{*}}$, которое для BPR-функций имеет вид $t \geqslant \bar {t}$, и лишь при $\mu \to 0{\kern 1pt} + $ минимум может выходить на границу ${{t}_{e}} = {{\bar {t}}_{e}}$)

(13)

${{\left\{ {{{\tau }_{e}}\left( {\bar {f}_{e}^{N}} \right)} \right\}}_{{e \in E}}} = \arg \mathop {\min }\limits_{t \in {\text{dom}}\sigma *} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \right\}.$

Из (14) следует, что

$\begin{gathered} \gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {t_{e}^{N}} \right)} \leqslant \\ \leqslant \;\mathop {\min }\limits_{t \in \operatorname{dom} \sigma *} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}{\kern 1pt} \left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi {\kern 1pt} \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \sum\limits_{e \in E} {\sigma _{e}^{*}{\kern 1pt} \left( {{{t}_{e}}} \right)} + \frac{1}{{2{{A}_{N}}}}\left\| {t - \bar {t}} \right\|_{2}^{2}} \right\} + \frac{\varepsilon }{2} \leqslant \\ \end{gathered} $

$\begin{gathered} \leqslant \mathop {\min }\limits_{t \in {\text{dom}}\sigma *} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \right\} + \\ + \;\frac{1}{{{{A}_{N}}}}\underbrace {\frac{1}{2}\sum\limits_{e \in E} {{{{\left( {{{\tau }_{e}}\left( {\bar {f}_{e}^{N}} \right) - {{{\bar {t}}}_{e}}} \right)}}^{2}}} }_{{{{\tilde {R}}}^{2}}} + \frac{\varepsilon }{2}. \\ \end{gathered} $

Следовательно,

$\begin{gathered} \gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {t_{e}^{N}} \right)} - \\ - \;\mathop {\min }\limits_{t \in {\text{dom}}\sigma *} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \right\} \leqslant \frac{{{{{\tilde {R}}}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}. \\ \end{gathered} $

Учитывая, что

$\begin{gathered} - \mathop {\min }\limits_{t \in \operatorname{dom} \sigma *} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t} \right\rangle } } \right] + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \right\} = \\ = \mathop {\max }\limits_{t \in \operatorname{dom} \sigma *} \left\{ {\left\langle {\frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}{{f}^{k}}} ,t} \right\rangle - \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {{{t}_{e}}} \right)} } \right\} = \sum\limits_{e \in E} {{{\sigma }_{e}}\left( {\bar {f}_{e}^{N}} \right)} , \\ \end{gathered} $

$\begin{gathered} - \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) - \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),{{y}^{k}}} \right\rangle } \right)} = \\ = \;\frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}} \left( {\left\langle {{{f}^{k}},{{y}^{k}}} \right\rangle + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{k}\ln \left( {{{x_{p}^{k}} \mathord{\left/ {\vphantom {{x_{p}^{k}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } - \left\langle {{{f}^{k}},{{y}^{k}}} \right\rangle } \right) = \\ = \;\gamma \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}\sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{k}\ln \left( {{{x_{p}^{k}} \mathord{\left/ {\vphantom {{x_{p}^{k}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \geqslant \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } , \\ \end{gathered} $

получаем следующую оценку сверху на зазор двойственности:

(14)

$\begin{gathered} 0 \leqslant \left\{ {\gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {t_{e}^{N}} \right)} } \right\} + \left\{ {\sum\limits_{e \in E} {{{\sigma }_{e}}\left( {\bar {f}_{e}^{N}} \right)} + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \right\} \leqslant \hfill \\ \leqslant \left\{ {\gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {t_{e}^{N}} \right)} } \right\} + \left\{ {\sum\limits_{e \in E} {{{\sigma }_{e}}\left( {\bar {f}_{e}^{N}} \right)} + \gamma \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}} \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{k}\ln \left( {{{x_{p}^{k}} \mathord{\left/ {\vphantom {{x_{p}^{k}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \right\} \leqslant \frac{{{{{\tilde {R}}}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}. \hfill \\ \end{gathered} $

Следовательно, имеем

(15)

$\begin{gathered} \sum\limits_{e \in E} {{{\sigma }_{e}}\left( {\bar {f}_{e}^{N}} \right)} + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } - \left( {\sum\limits_{e \in E} {{{\sigma }_{e}}(f_{e}^{*})} + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{*}\ln \left( {{{x_{p}^{*}} \mathord{\left/ {\vphantom {{x_{p}^{*}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \right) \leqslant \\ \leqslant \;\left\{ {\gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \sum\limits_{e \in E} {\sigma _{e}^{*}\left( {t_{e}^{N}} \right)} } \right\} + \left\{ {\sum\limits_{e \in E} {{{\sigma }_{e}}\left( {\bar {f}_{e}^{N}} \right)} + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \right\} \leqslant \frac{{{{{\tilde {R}}}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}, \\ \end{gathered} $

где $\left( {f*,x{\text{*}}} \right)$ – решение задачи (1).

На формулу (14) можно смотреть как на критерий останова УМПТ. А именно, ждем когда (вычисляемый с помощью быстрого автоматического дифференцирования п. 6) зазор двойственности станет меньше $\varepsilon $. Формула (10) (с заменой $R$ на $\tilde {R}$) содержит гарантированную оценку, на число итераций УМПТ, после которого метод гарантированно остановится по критерию (14).

Для модели (стохастической) стабильной динамики ($\mu \to 0{\kern 1pt} + $) необходимо немного по-другому провести рассуждения:

$\begin{gathered} \gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\bar {f},{{t}^{N}} - \bar {t}} \right\rangle \leqslant \\ \leqslant \;\mathop {\min }\limits_{t \geqslant \bar {t}} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \left\langle {\bar {f},t - \bar {t}} \right\rangle + \frac{1}{{2{{A}_{N}}}}\left\| {t - \bar {t}} \right\|_{2}^{2}} \right\} + \frac{\varepsilon }{2} \leqslant \\ \leqslant \;\mathop {\min }\limits_{t \geqslant \bar {t},\;\left\| {t - \bar {t}} \right\|_{2}^{2} \leqslant 2{{R}^{2}}} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left( {\gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t - {{y}^{k}}} \right\rangle } \right)} } \right] + \left\langle {\bar {f},t - \bar {t}} \right\rangle } \right\} + \frac{{{{R}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}, \\ \end{gathered} $

где ${{R}^{2}} = \frac{1}{2}\left\| {{{t}_{*}} - {{y}^{0}}} \right\|_{2}^{2} = \frac{1}{2}\left\| {{{t}_{*}} - \bar {t}} \right\|_{2}^{2}$. Следовательно, имеем

Поскольку

$\begin{gathered} - \mathop {\min }\limits_{t \geqslant \bar {t},\;\left\| {t - \bar {t}} \right\|_{2}^{2} \leqslant 10{{R}^{2}}} \left\{ {\frac{1}{{{{A}_{N}}}}\left[ {\sum\limits_{k = 0}^N {{{\alpha }_{k}}\left\langle {\nabla \gamma \psi \left( {{{{{y}^{k}}} \mathord{\left/ {\vphantom {{{{y}^{k}}} \gamma }} \right. \kern-0em} \gamma }} \right),t} \right\rangle } } \right] + \left\langle {\bar {f},t - \bar {t}} \right\rangle } \right\} = \\ = \;\mathop {\max }\limits_{t \geqslant \bar {t},\;\left\| {t - \bar {t}} \right\|_{2}^{2} \leqslant 10{{R}^{2}}} \left\{ {\left\langle {\frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}{{f}^{k}}} ,t} \right\rangle - \left\langle {\bar {f},t - \bar {t}} \right\rangle } \right\} = \\ = \;\mathop {\max }\limits_{t \geqslant \bar {t},\;\left\| {t - \bar {t}} \right\|_{2}^{2} \leqslant 10{{R}^{2}}} \left\{ {\left\langle {{{{\bar {f}}}^{N}} - \bar {f},t - \bar {t}} \right\rangle } \right\} + \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle \geqslant \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle + 3R{{\left\| {{{{\left( {{{{\bar {f}}}^{N}} - \bar {f}} \right)}}_{ + }}} \right\|}_{2}}, \\ \end{gathered} $

то

(16)

$\begin{gathered} \gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\bar {f},{{t}^{N}} - \bar {t}} \right\rangle + \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } + 3R{{\left\| {{{{\left( {{{{\bar {f}}}^{N}} - \bar {f}} \right)}}_{ + }}} \right\|}_{2}} \leqslant \\ \leqslant \gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\bar {f},{{t}^{N}} - \bar {t}} \right\rangle + \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle + \gamma \frac{1}{{{{A}_{N}}}}\sum\limits_{k = 0}^N {{{\alpha }_{k}}} \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{N}\ln \left( {{{x_{p}^{N}} \mathord{\left/ {\vphantom {{x_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } + 3R{{\left\| {{{{\left( {{{{\bar {f}}}^{N}} - \bar {f}} \right)}}_{ + }}} \right\|}_{2}} \leqslant \frac{{5{{R}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}. \\ \end{gathered} $

Повторяя рассуждения разд. 3 (см. [33]) и разд. 6.11 (см. [36]) из (16), получаем

(17)

$\begin{gathered} \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } - \left( {\left\langle {f*,\bar {t}} \right\rangle + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {x_{p}^{*}\ln \left( {{{x_{p}^{*}} \mathord{\left/ {\vphantom {{x_{p}^{*}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } } \right) \leqslant \\ \leqslant \;\gamma \psi \left( {{{{{t}^{N}}} \mathord{\left/ {\vphantom {{{{t}^{N}}} \gamma }} \right. \kern-0em} \gamma }} \right) + \left\langle {\bar {f},{{t}^{N}} - \bar {t}} \right\rangle + \left\langle {{{{\bar {f}}}^{N}},\bar {t}} \right\rangle + \gamma \sum\limits_{w \in OD} {\sum\limits_{p \in {{P}_{w}}} {\bar {x}_{p}^{N}\ln \left( {{{\bar {x}_{p}^{N}} \mathord{\left/ {\vphantom {{\bar {x}_{p}^{N}} {{{d}_{w}}}}} \right. \kern-0em} {{{d}_{w}}}}} \right)} } \leqslant \frac{{5{{R}^{2}}}}{{{{A}_{N}}}} + \frac{\varepsilon }{2}, \\ \end{gathered} $

(18)

${{\left\| {{{{\left( {{{{\bar {f}}}^{N}} - \bar {f}} \right)}}_{ + }}} \right\|}_{2}} \leqslant \frac{{5R}}{{{{A}_{N}}}} + \frac{\varepsilon }{{2R}},$

где $\left( {f*,x{\text{*}}} \right)$ – решение задачи (2).

На формулу (16) можно смотреть как на критерий останова УМПТ. А именно, ждем когда (вычисляемая) левая часть (16) станет меньше $\varepsilon $. Формула (10) (с заменой $R$ на $\sqrt 5 R$) содержит гарантированную оценку, на число итераций УМПТ, после которого метод гарантированно остановится по критерию (16).

Далее для краткости будем единым образом обозначать через $\bar {R}$ либо $\tilde {R}$, либо $\sqrt 5 R$. Расшифровывать $\bar {R}$ нужно будет в зависимости от контекста. Если $\mu \to 0{\kern 1pt} + $, то $\bar {R} = \sqrt 5 R$, иначе $\bar {R} = \tilde {R}$.

Все приведенные здесь формулы выдерживают предельные переходы $\gamma \to 0{\kern 1pt} + $, $\mu \to 0{\kern 1pt} + $ (см. разд. 2).

6. ВЫЧИСЛЕНИЕ (СУБ-)ГРАДИЕНТОВ В ЗАДАЧЕ ПОИСКА РАВНОВЕСИЙ В ТРАНСПОРТНЫХ СЕТЯХ

Приведем, следуя [5], [27], сглаженный идемпотентный аналог метода Форда–Беллмана (см. [29], [41]), позволяющий эффективно рассчитывать значение характеристической функции $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$. Для этого предположим, что любые движения по ребрам графа с учетом их ориентации являются допустимыми, т.е. множество путей, соединяющих заданные две вершины (источник и сток), – это множество всевозможных способов добраться из источника в сток по ребрам имеющегося графа с учетом их ориентации. Этого всегда можно добиться раздутием исходного графа в несколько раз за счет введения дополнительных вершин и ребер. Такое раздутие заведомо можно сделать за ${\rm O}\left( n \right)$. Отметим, что при этом в качестве путей будут присутствовать, в том числе, и самопересекающиеся маршруты. Однако можно показать, что вклад таких “не физических” путей в итоговую равновесную конфигурацию будет пренебрежимо мал.

Будем считать, как и прежде, что число ребер в любом пути не больше $H = {\rm O}\left( {\sqrt n } \right)$. Введем классы путей: $P_{{ij}}^{l}$ – множество всех путей из $i$ в $j$, состоящих ровно из $l$ ребер, $\tilde {P}_{{ij}}^{l}$ – множество всех путей из $i$ в $j$, состоящих из не более чем $l$ ребер. Зафиксируем источник (вершину) $i \in V$, и введем следующие функции для $j \in V$, $l = 1,\;...,\;H$:

$a_{{ij}}^{l}\left( t \right) = \gamma {{\psi }_{{P_{{ij}}^{l}}}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = \gamma \ln \left( {\sum\limits_{p \in P_{{ij}}^{l}} {\exp \left( { - {{\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \mathord{\left/ {\vphantom {{\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \gamma }} \right. \kern-0em} \gamma }} \right)} } \right),$

$b_{{ij}}^{l}\left( t \right) = \gamma {{\psi }_{{\tilde {P}_{{ij}}^{l}}}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = \gamma \ln \left( {\sum\limits_{p \in \tilde {P}_{{ij}}^{l}} {\exp \left( { - {{\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \mathord{\left/ {\vphantom {{\sum\limits_{e \in E} {{{\delta }_{{ep}}}{{t}_{e}}} } \gamma }} \right. \kern-0em} \gamma }} \right)} } \right).$

Некоторые из этих функций могут быть равны $ - \infty $. Это означает, что соответствующее множество маршрутов пустое. Данные функции можно вычислять рекурсивным образом:

(19)

$\begin{gathered} a_{{ij}}^{1}\left( t \right) = b_{{ij}}^{1}\left( t \right) = \left\{ \begin{gathered} - {{t}_{e}},\quad e = \left( {i \to j} \right) \in E, \hfill \\ - \infty \quad e = \left( {i \to j} \right) \notin E, \hfill \\ \end{gathered} \right. \\ a_{{ij}}^{{l + 1}}\left( t \right) = \gamma \ln \left( {\sum\limits_{k:\;e = \left( {k \to j} \right) \in E} {\exp \left( {{{\left( {a_{{ik}}^{l}\left( t \right) - {{t}_{e}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{ik}}^{l}\left( t \right) - {{t}_{e}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)} } \right), \\ b_{{ij}}^{{l + 1}}\left( t \right) = \gamma \ln \left( {\exp \left( {{{b_{{ij}}^{l}\left( t \right)} \mathord{\left/ {\vphantom {{b_{{ij}}^{l}\left( t \right)} \gamma }} \right. \kern-0em} \gamma }} \right) + \exp \left( {{{a_{{ij}}^{{l + 1}}\left( t \right)} \mathord{\left/ {\vphantom {{a_{{ij}}^{{l + 1}}\left( t \right)} \gamma }} \right. \kern-0em} \gamma }} \right)} \right),\quad j \in V,\quad l = 1,\;...,\;H - 1. \\ \end{gathered} $

На каждом шаге $l$ необходимо сделать ${\rm O}\left( n \right)$ арифметических операций. Следовательно, для вычисления $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ необходимо сделать ${\rm O}\left( {SHn} \right)$ арифметических операций, где $S = \left| O \right|$ – число источников, как правило, можно считать $S \ll n$. Причем вычисление функции $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ (и ее градиента) может быть распараллелено на $S$ процессорах. При $\gamma \to 0{\kern 1pt} + $ процедура вырождается в известный метод Форда–Беллмана (см. [29], [41]) (динамическое программирование). В процедуре Форда–Беллмана требуется посчитать $H$-степень матрицы $A = {{\left\| {{{a}_{{ij}}}} \right\|}_{{i,j \in V}}}$,

${{a}_{{ij}}} = {{t}_{e}},\quad e = \left( {i \to j} \right) \in E,$

${{a}_{{ij}}} = \infty ,\quad e = \left( {i \to j} \right) \notin E,$

в идемпотентной математике (вместо обычного поля используется тропическое полуполе [42] со следующими операциями: сложение $a \oplus b = \min \left\{ {a,b} \right\}$, произведение $a \otimes b = a + b$). Учитывая, что

$a \oplus b = \min \left\{ {a,b} \right\} = - \mathop {\lim }\limits_{\gamma \to 0 + } \gamma \ln \left( {\exp \left( { - {a \mathord{\left/ {\vphantom {a \gamma }} \right. \kern-0em} \gamma }} \right) + \exp \left( { - {b \mathord{\left/ {\vphantom {b \gamma }} \right. \kern-0em} \gamma }} \right)} \right),$

$a \otimes b = a + b = - \mathop {\lim }\limits_{\gamma \to 0 + } \gamma \ln \left( {\exp \left( { - {a \mathord{\left/ {\vphantom {a \gamma }} \right. \kern-0em} \gamma }} \right)\exp \left( { - {b \mathord{\left/ {\vphantom {b \gamma }} \right. \kern-0em} \gamma }} \right)} \right),$

можно посчитать обычную (над обычным полем) $H$-степень матрицы ${{A}^{\gamma }} = {{\left\| {a_{{ij}}^{\gamma }} \right\|}_{{i,j \in V}}}$,

$a_{{ij}}^{\gamma } = {{e}^{{ - {{{{t}_{e}}} \mathord{\left/ {\vphantom {{{{t}_{e}}} \gamma }} \right. \kern-0em} \gamma }}}},\quad e = \left( {i \to j} \right) \in E,$

$a_{{ij}}^{\gamma } = 0,\quad e = \left( {i \to j} \right) \notin E,$

и применить поэлементно к полученной матрице $ - \gamma \ln \left( \cdot \right)$. В пределе $\gamma \to 0{\kern 1pt} + $ получим метод Форда–Беллмана. Однако если не делать предельный переход, то получается нужный нам сглаженный вариант этого алгоритма с такой же временной сложностью. Некоторый аналог этого сглаженного варианта, по сути, и был описан выше.

Используя быстрое автоматическое дифференцирование (см. [28]), опишем способ вычисления градиента функции $\gamma {{\psi }^{i}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = \sum\nolimits_{j:\;w = \left( {i,j} \right) \in OD} {{{d}_{w}}\gamma {{\psi }_{w}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)} $:

(20)

$\begin{gathered} \frac{{\partial {{\psi }^{i}}}}{{\partial b_{{ij}}^{H}}} = {{d}_{{ij}}},\quad \frac{{\partial {{\psi }^{i}}}}{{\partial a_{{ij}}^{H}}} = 0,\quad \frac{{\partial {{\psi }^{i}}}}{{\partial b_{{ij}}^{l}}} = \frac{{\partial {{\psi }^{i}}}}{{\partial b_{{ij}}^{{l + 1}}}}\frac{{\partial b_{{ij}}^{{l + 1}}}}{{\partial b_{{ij}}^{l}}}, \\ \frac{{\partial {{\psi }^{i}}}}{{\partial a_{{ij}}^{l}}} = \frac{{\partial {{\psi }^{i}}}}{{\partial b_{{ij}}^{l}}}\frac{{\partial b_{{ij}}^{l}}}{{\partial a_{{ij}}^{l}}} + \sum\limits_{k:\;e = \left( {k \to j} \right) \in E} {\frac{{\partial {{\psi }^{i}}}}{{\partial a_{{ik}}^{{l + 1}}}}\frac{{\partial a_{{ik}}^{{l + 1}}}}{{\partial a_{{ik}}^{l}}}} ,\quad j \in V,\quad l = H - 1,\;...,\;1. \\ \end{gathered} $

(21)

$\frac{{\partial {{\psi }^{i}}}}{{\partial {{t}_{e}}}} = \sum\limits_{l = 0}^{H - 1} {\sum\limits_{j \in V} {\frac{{\partial {{\psi }^{i}}}}{{\partial a_{{ij}}^{{l + 1}}}}\frac{{\partial a_{{ij}}^{{l + 1}}}}{{\partial {{t}_{e}}}} = } } \sum\limits_{l = 0}^{H - 1} {\frac{{\partial {{\psi }^{i}}}}{{\partial a_{{ij'}}^{{l + 1}}}}\frac{{\partial a_{{ij'}}^{{l + 1}}}}{{\partial {{t}_{e}}}},\quad e = \left( {k,\;j'} \right)} .$

Частные производные

$\frac{{\partial b_{{ij}}^{{l + 1}}}}{{\partial b_{{ij}}^{l}}},\quad \frac{{\partial b_{{ij}}^{l}}}{{\partial a_{{ij}}^{l}}},\quad \frac{{\partial a_{{ik}}^{{l + 1}}}}{{\partial a_{{ik}}^{l}}},\quad \frac{{\partial a_{{ij}}^{{l + 1}}}}{{\partial {{t}_{e}}}}$

могут быть явно вычислены из системы (19) за ${\rm O}\left( 1 \right)$ каждая. Поясним примером. Пусть, например, для некоторых $l$ и $j$ имеет место

$a_{{ij}}^{{l + 1}}\left( t \right) = \gamma \ln \left( {\exp \left( {{{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right) + \exp \left( {{{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)} \right).$

Тогда верно следующее:

$\frac{{\partial a_{{ij}}^{{l + 1}}}}{{\partial a_{{i{{k}_{1}}}}^{l}}} = \frac{{\exp \left( {{{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)}}{{\exp \left( {{{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right) + \exp \left( {{{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)}},$

$\frac{{\partial a_{{ij}}^{{l + 1}}}}{{\partial {{t}_{{{{e}_{1}}}}}}} = - \frac{{\exp \left( {{{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)}}{{\exp \left( {{{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{1}}}}^{l}\left( t \right) - {{t}_{{{{e}_{1}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right) + \exp \left( {{{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \mathord{\left/ {\vphantom {{\left( {a_{{i{{k}_{2}}}}^{l}\left( t \right) - {{t}_{{{{e}_{2}}}}}} \right)} \gamma }} \right. \kern-0em} \gamma }} \right)}}.$

(Минимум в (13) достигается во внутренней точке множества $\operatorname{dom} \sigma {\text{*}}$, которое для BPR-функций имеет вид $t \geqslant \bar {t}$, и лишь при $\mu \to 0{\kern 1pt} + $ минимум может выходить на границу ${{t}_{e}} = {{\bar {t}}_{e}}$.)

Остальные частные производные являются неизвестными. Несложно понять, что система (20) может быть последовательно разрешена за ${\rm O}\left( {Hn} \right)$. Фактически тут происходит процесс прохождения графа вычислений (19) в обратном направлении, с той лишь разницей, что обратный проход получается приблизительно в 4 раза дороже. К сожалению, в отличие от прямого процесса (19), теперь, чтобы вычислить градиент, необходимо хранить в памяти промежуточные вычисления, см. формулу (21),т.е. для вычисления градиента $\gamma {{\psi }^{i}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ требуются время ${\rm O}\left( {Hn} \right)$ и память ${\rm O}\left( {Hn} \right)$. Таким образом, для вычисления градиента $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right) = \sum\nolimits_{i \in O} {\gamma {{\psi }^{i}}\left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)} $ требуется время ${\rm O}\left( {SHn} \right)$ и память ${\rm O}\left( {SHn} \right)$. Также как и при вычислении значения функции $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$, описанная выше схема вычисления градиента может быть распараллелена на $S = \left| O \right|$ процессорах.

Заметим, что расчет $\nabla \gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$, реализованный на языке программирования MatLab 8 на стандартном современном ноутбуке (с тактовой частотой 1.9 ГГц), для манхетенского графа c числом ребер $n \simeq {{10}^{3}}$, занимал порядка одной минуты, в то время как использование прямого дифференцирования занимало порядка 100 мин.

В связи с написанным выше, важно отметить, в пределе $\gamma \to 0{\kern 1pt} + $ оценка ${\rm O}\left( {SHn} \right)$ сложности вычисления градиента $\gamma \psi \left( {{t \mathord{\left/ {\vphantom {t \gamma }} \right. \kern-0em} \gamma }} \right)$ может быть улучшена до оценки ${\rm O}\left( {Sn\ln n} \right)$ вычисления субградиента (5), (6) [6], [29]. Действительно, для каждого из $S$ источников можно построить (например, алгоритмом Дейкстры, см. [29]) соответствующее дерево кратчайших путей. Исходя из принципа динамического программирования “часть кратчайшего пути сама будет кратчайшим путем” несложно понять, что получится именно дерево, с корнем в рассматриваемом источнике. Это можно сделать для одного источника за ${\rm O}\left( {n\ln n} \right)$ [6], [12], [29]. Однако, главное, правильно взвешивать ребра (их $n$) такого дерева, чтобы за один проход этого дерева можно было восстановить вклад (по всем ребрам) соответствующего источника в субградиент. Ребро должно иметь вес, равный сумме всех проходящих через него корреспонденций с заданным источником (корнем дерева). Имея значения соответствующих корреспонденций (их не больше ${\rm O}\left( n \right)$) за один обратный проход (т.е. с листьев к корню) такого дерева можно осуществить необходимое взвешивание (с затратами не более ${\rm O}\left( n \right)$). Делается это по правилу: вес ребра равен сумме корреспонденции (возможно, равной нулю), в соответствующую вершину, в которую ребро входит, и сумме весов всех ребер (если таковые имеются), выходящих из упомянутой вершины.

7. СОПОСТАВЛЕНИЕ ОЦЕНОК, ЧИСЛЕННЫЕ ЭКСПЕРИМЕНТЫ

Используя наработки предыдущих разделов, в частности, формулы (10)–(12), (15), (17), (18), можно резюмировать полученные в статье результаты в виде следующей теоремы.

Теорема. УМПТ из разд. 4 с критериями останова (14), (16) гарантированно остановится, достигнув желаемой точности $\varepsilon $, сделав арифметических операций не больше, чем приведено в соответствующем поле таблицы 1.

Таблица 1.

Время работы	$\gamma > 0$	$\gamma \to 0 + $
$\mu \to 0{\kern 1pt} + $, $\mu > 0$ – не важно	${\rm O}\left( {SHn\sqrt {\frac{{Hd{{{\bar {R}}}^{2}}}}{{\gamma \varepsilon }}} } \right)$ (22)	${\rm O}\left( {Sn\ln n\frac{{H{{d}^{2}}{{{\bar {R}}}^{2}}}}{{{{\varepsilon }^{2}}}}} \right)$ (23)

Обратим внимание, что оценки (22), (23) не зависят от (потенциально экспоненциально большого) числа путей $\left| P \right|$. В частности, $\left| P \right| \gg {{2}^{{\sqrt n }}}$, для манхетенских сетей. Также обратим внимание, что оценка (22) весьма чувствительна к предельному переходу $\gamma \to 0 + $.

Важно также заметить, что обычные (не стохастические) равновесия можно искать с помощью искусственного введения энтропийной регуляризации (см. [9], [17]). При этом, чтобы с точностью $\varepsilon > 0$ по функции решить исходную задачу (1) или (2) с $\gamma = 0$, можно действовать следующим образом. Выбрать

${{\gamma }_{*}} = \frac{\varepsilon }{{2\sum\limits_{w \in OD} {{{d}_{w}}\ln \left| {{{P}_{w}}} \right|} }}$

и решать регуляризованную задачу (1) или (2) с $\gamma = {{\gamma }_{*}}$, но с точностью ${\varepsilon \mathord{\left/ {\vphantom {\varepsilon 2}} \right. \kern-0em} 2}$ (см. [17]). В этом случае оценка (22) из таблицы 1 примет вид

(24)

${\rm O}\left( {SHn\frac{d}{\varepsilon }\sqrt {H{{{\bar {R}}}^{2}}\omega } } \right),$

где $\omega = \mathop {\max }\limits_{w \in OD} \ln \left| {{{P}_{w}}} \right|$ (в худшем случае можно ожидать $\omega \sim \sqrt n $). Сопоставляя оценки (23), (24), можно прийти к выводу, что введенная регуляризация (приводящая к сглаживанию двойственной задачи (3), см. [38]) оправдана. Однако не стоит забывать, что все приведенные оценки – это верхние оценки. Численные эксперименты показали, что эти оценки не всегда точны, особенно в части формулы (23).

Прежде всего, заметим, что оценки общего числа арифметических операций из таблицы 1 можно понимать как [стоимость итерации] $ \times $ [число итераций]. А число итераций $N$, с некоторой натяжкой, можно просто считать пропорциональным относительной точности $\tilde {\varepsilon }$ в некоторой степени $N \sim {{\tilde {\varepsilon }}^{{ - \beta }}}$. ($\tilde {\varepsilon } = 1\% = 0.01$ – означает, что начальную невязку по функции или зазору двойственности необходимо уменьшить в 100 раз.)

Именно такая (относительная) точность, как правило, интересна на практике. Из теоретических верхних оценок (22)–(24) напрашивается вывод, что при $\gamma \gg {{\gamma }_{*}}$ имеем $N \sim {{\tilde {\varepsilon }}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0em} 2}}}}$, при $\gamma \simeq {{\gamma }_{*}}$ имеем $N \sim {{\tilde {\varepsilon }}^{{ - 1}}}$, а при $\gamma \to 0{\kern 1pt} + $ имеем $N \sim {{\tilde {\varepsilon }}^{{ - 2}}}$. На самом деле, конечно, опущенные числовые множители тут также играют важную роль. Однако еще более важно то, что численные эксперименты (см. [43]) не подтвердили оценку $N \sim {{\tilde {\varepsilon }}^{{ - 2}}}$ при $\gamma \to 0 + $. Более того, наблюдалась совсем другая зависимость $N \sim {{\tilde {\varepsilon }}^{{ - 1}}}$ (см. [43]), точнее говоря, наблюдалась зависимость $N \sim {{С }_{1}} + {{С }_{2}}{{\tilde {\varepsilon }}^{{ - 1}}}$ (см. фиг. 1, 2 ). Учитывая, что стоимость итерации заметно меньше в случае $\gamma \to 0{\kern 1pt} + $ ($SHn \to Sn\ln n$, см. разд. 6), то вывод об оправданности регуляризации приходится подставить под сомнение. Для графа города Анахайм (Anaheim) [44] с $n \simeq {{10}^{3}}$, $S \sim 40$ и $\left| {OD} \right| \sim 1.5 \times {{10}^{3}}$ УМПТ, реализованный на стандартном современном ноутбуке (с тактовой частотой 1.9 ГГц) на языке программирования Python 2.4 [43], находил с относительной точностью $\tilde {\varepsilon } \simeq 0.01$ равновесие в модели Бэкмана ($\mu > 0$, $\gamma \to 0{\kern 1pt} + $) приблизительно за 7 мин, в модели стабильной динамики ($\mu > 0$, $\gamma \to 0{\kern 1pt} + $) за 10 мин (см. фиг. 1, 2 ). Для соответствующей регуляризованной модели – заметно дольше [45].

Фиг. 1.

Модель Бэкмана ($\mu = 0.25$), $\gamma \to 0{\kern 1pt} + $ ($\tilde {\varepsilon }$ обозначено здесь через $\varepsilon $).

Фиг. 2.

Модель стабильной динамики ($\mu \to 0{\kern 1pt} + $), $\gamma \to 0{\kern 1pt} + $ ($\tilde {\varepsilon }$ обозначено здесь через $\varepsilon $).

Отметим, что написанный код не был оптимизирован. Мы полагаем, что при правильной реализации и с использованием языка более низкого уровня (например, С++) можно ожидать выигрыша как минимум на порядок (скорее всего, на два). В таком случае можно рассчитывать, что предлагаемый в данной статье метод приблизится по производительности к методу условного градиента [3], [6] (Франк–Вульфа), который считается сейчас наиболее эффективным методом поиска равновесного распределения потоков по путям (ребрам) в модели Бэкмана ($\mu = 0.25$, $\gamma \to 0 + $). В лучшей (из известных нам) реализаций этого метода, сделанной А.С. Аникиным и А.Ю. Горновым на С++ (см. [45]), для графа, с аналогичными параметрами и с аналогичными требованиями к точности, метод условного градиента сходился за пару секунд. Для модели стабильной динамики имеющиеся сейчас реализации (при $\gamma \to 0{\kern 1pt} + $) существующих методов (см. [6], [12]) заметно проигрывают изложенным в данной работе.

Описанный в статье подход можно переписать, беря в качестве базисного метода УМПТ из статьи [17], в варианте, работающем с сильно выпуклыми постановками задач. Поскольку двойственная задача (3) не является сильно выпуклой (во всяком случае, на данный момент умеют устанавливать только выпуклость этой задачи), то необходимо регуляризовывать двойственную задачу. Можно показать, по сути, рассуждая подобно статье [39] (см. также концовку статьи [33]), что при правильной регуляризации существенно упрощаются формулы восстановления решения прямой задачи.

Описанный в статье формализм может быть перенесен и на задачи поиска равновесий в многостадийных транспортных моделях [1], [7], [8], [10], [11], [13], [30]. Для этого стоит использовать подход из работ [11], [13].

Авторы выражают благодарность Ю.Е. Нестерову за ряд ценных замечаний.

Список литературы

Гасников А.В., Дорн Ю.В., Нестеров Ю.Е., Шпирко С.В. О трехстадийной версии модели стационарной динамики транспортных потоков // Матем. моделирование. 2014. Т. 26. № 6. С. 34–70. arXiv:1405.7630
Sheffi Y. Urban transportation networks: Equilibrium analysis with mathematical programming methods. N.J.: Prentice–Hall Inc., Englewood Cliffs, 1985.
Patriksson M. The traffic assignment problem. Models and methods. Utrecht, Netherlands: VSP, 1994.
Nesterov Yu., de Palma A. Stationary dynamic solutions in congested transportation Networks: Summary and Perspectives // Networks Spatial Econ. 2003. № 3(3). P. 371–395.
Введение в математическое моделирование транспортных потоков / Под ред. А.В. Гасникова. Изд. 2. М.: МЦНМО, 2013. 427 с. http://www.mou.mipt.ru/gasnikov1129.pdf
Гасников А.В., Двуреченский П.Е., Дорн Ю.В., Максимов Ю.В. Численные методы поиска равновесного распределения потоков в модели Бэкмана и модели стабильной динамики // Матем. моделирование. 2016. Т. 28. № 10. С. 40–64. arXiv:1506.00293
Бабичева Т.С., Гасников А.В., Лагуновская А.А., Мендель М.А. Двухстадийная модель равновесного распределения транспортных потоков // Труды МФТИ. 2015. Т. 7. № 3. С. 31–41. https://mipt.ru/upload/medialibrary/971/31-41.pdf
Гасников А.В. Об эффективной вычислимости конкурентных равновесий в транспортно-экономических моделях // Матем. моделирование. 2015. Т. 27. № 12. С. 121–136. arXiv:1410.3123
Гасников А.В., Гасникова Е.В., Двуреченский П.Е., Ершов Е.И., Лагуновская А.А. Поиск стохастических равновесий в транспортных моделях равновесного распределения потоков // Труды МФТИ. 2015. Т. 7. № 4. С. 114–128. https://mipt.ru/upload/medialibrary/93f/114-128.pdf
Гасников А.В., Двуреченский П.Е., Камзолов Д.И., Нестеров Ю.Е., Спокойный В.Г., Стецюк П.И., Суворикова А.Л., Чернов А.В. Поиск равновесий в многостадийныйх транспортных моделях // Труды МФТИ. 2015. Т. 7. № 4. С. 143–155. https://mipt.ru/upload/medialibrary/ffe/143-155.pdf
Гасников А.В., Гасникова Е.В., Мациевский С.В., Усик И.В. О связи моделей дискретного выбора с разномасштабными по времени популяционными играми загрузок // Труды МФТИ. 2015. Т. 7. № 4. С. 129–142.
Аникин А.С., Гасников А.В., Горнов А.Ю., Двуреченский П.Е., Семенов В.В. Параллелизуемые двойственные методы поиска равновесий в смешанных моделях распределения потоков в больших транспортных сетях. Труды 40-й междунар. школы-конф. “Информационные технологии и системы–2016”. Россия, С.-Петербург (Репино), 25–30 сентября 2016. 8 с. arXiv:1604.08183
Dvurechensky P., Gasnikov A., Gasnikova E., Matsievski S., Rodomanov A., Usik I. Primal-Dual Method for Searching Equilibrium in Hierarchical Congestion Population Games. Supplementary Proc. of the 9th Internat. Conference on Discrete Optimizat. and Operat. Research and Sci. School (DOOR-2016). Vladivostok, 19–23 September, 2016. P. 584–595. arXiv:1606.08988
Гасников А.В. Эффективные численные методы поиска равновесий в больших транспортных сетях. Дис. д. физ.-матем. наук. 05.13.18. М.: МФТИ, 2016. 487 с. arXiv:1701.00719
Nesterov Yu. Universal gradient methods for convex optimization problems // CORE Discussion Paper 2013/63. 2013. http://www.uclouvain.be/cps/ucl/doc/core/documents/coredp2013_26web.pdf
Nesterov Yu. Universal gradient methods for convex optimization problems // Math. Program. Ser. A. 2015. V. 152. P. 381–404.
Гасников А.В., Нестеров Ю.Е. Универсальный метод для задач стохастической композитной оптимизации // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 7. С. в печати arXiv:1604.05275
Немировский А.С., Юдин Д.Б. Сложность задач и эффективность методов оптимизации. М.: Наука, 1979. 384 с.
Sandholm W. Population games and Evolutionary dynamics. Economic Learning and Social Evolution. MIT Press; Cambridge, 2010.
Andersen S.P., de Palma A., Thisse J.-F. Discrete choice theory of product differentiation. MIT Press; Cambridge, 1992.
Ethier N.S., Kurtz T.G. Markov processes. Wiley Series in Probability and Mathematical Statistics: Probability and Mathematical Statistics. New York: Wiley, 1986.
Малышев В.А., Пирогов С.А. Обратимость и необратимость в стохастической химической кинетике // Успехи матем. наук. 2008. Т. 63. Вып. 1(379). С. 4–36.
Гасников А.В., Гасникова Е.В. Об энтропийно-подобных функционалах, возникающих в стохастической химической кинетике при концентрации инвариантной меры и в качестве функций Ляпунова динамики квазисредних // Матем. заметки. 2013. Т. 94. № 6. С. 816–824. arXiv:1410.3126
Баймурзина Д.Р., Гасников А.В., Гасникова Е.В. Теория макросистем с точки зрения стохастической химической кинетики // Труды МФТИ. 2015. Т. 7. № 4. С. 95–103. https://mipt.ru/upload/medialibrary/ae2/95-103.pdf
Санов И.Н. О вероятности больших отклонений случайных величин // Матем. сб. 1957. Т. 42(84). № 1. С. 11–44.
Бузун Н.О., Гасников А.В., Гончаров Ф.О. Горбачев О.Г., Гуз С.А., Крымова Е.А., Натан А.А., Черноусова Е.О. Стохастический анализ в задачах. Учебное пособие. Ч. 1 / Под ред. А.В. Гасникова. М.: МФТИ, 2016. 212 с. arXiv:1508.03461
Nesterov Y. Characteristic functions of directed graphs and applications to stochastic equilibrium problems // Optim. Engng. 2007. V. 8. P. 193–214.
Ким К., Нестеров Ю., Скоков В., Черкасский Б. Эффективные алгоритмы для дифференцирования и задачи экстремали // Экономика и матем. методы. 1984. Т. 20. С. 309–318.
Ahuja R.K., Magnati T.L., Orlin J.B. Network flows: Theory, algorithms and applications. Prentice Hall, 1993.
Ortúzar J.D., Willumsen L.G. Modelling transport. New York: Willey, 2011.
Boyd S., Vandenberghe L. Convex optimization. Cambridge University Press, 2004.
Лидбеттер М., Линдгрен Г., Ротсен Х. Экстремумы случайных последовательностей и процессов. М.: Мир, 1989.
Аникин А.С., Гасников А.В., Двуреченский П.Е., Тюрин А.И., Чернов А.В. Двойственные подходы к задачам минимизации сильно выпуклых функционалов простой структуры при аффинных ограничениях // Ж. вычисл. матем. и матем. физ. 2017. Т. 57. № 6. arXiv:1602.01686
Chernov A., Dvurechensky P., Gasnikov A. Fast Primal-Dual Gradient Method for Strongly Convex Minimization Problems with Linear Constraints // In: Kochetov, Yu. et all (eds.) DOOR-2016. LNCS, V. 9869. P. 391–403. Heidelberg: Springer, 2016. arXiv:1605.02970
Nesterov Yu. Primal-dual subgradient methods for convex problems // Math. Program. Ser. B. 2009. V. 120(1). P. 261–283.
Nemirovski A., Onn S., Rothblum U.G. Accuracy certificates for computational problems with convex structure // Math. of Operat. Research. 2010. V. 35. № 1. P. 52–78.
Nesterov Yu. Gradient methods for minimizing composite functions // Math. Prog. 2013. V. 140. № 1. P. 125–161. http://www.uclouvain.be/cps/ucl/doc/core/documents/Composit.pdf
Nesterov Yu. Smooth minimization of non-smooth function // Math. Program. Ser. A. 2005. V. 103. № 1. P. 127–152.
Гасников А.В., Гасникова Е.В., Нестеров Ю.Е., Чернов А.В. Об эффективных численных методах решения задач энтропийно-линейного программирования // Ж. вычисл. матем. и матем. физ. 2016. Т. 56. № 4. С. 523–534. arXiv:1410.7719
Anikin A., Dvurechensky P., Gasnikov A., Golov A., Gornov A., Maximov Yu., Mendel M., Spokoiny V. Modern efficient numerical approaches to regularized regression problems in application to traffic demands matrix calculation from link loads // Proc. of Internat.Conference ITAS-2015. Russia, Sochi, September, 2015. arXiv:1508.00858
Форд Л., Фалкерсон Д. Потоки в сетях. М.: Мир, 1966.
Kolokoltsov V.N., Maslov V.P. Idempotent analysis and applications. Dordrecht: Kluwer, 1997.
https://github.com/dilyararimovna/TransportFlows2016
https://github.com/bstabler/TransportationNetworks
Аникин А.С. Программная реализация метода условного градиента Франк–Вульфа для поиска равновесия в модели Бэкмана. Выступление на отчетной транспортной конференции 26 декабря 2015 г. МЦНМО. http://www.mathnet.ru/php/presentation.phtml?option_lang=rus&presentid=13273

Дополнительные материалы отсутствуют.

Инструменты

следующая статья выпуска предыдущая статья выпуска содержание выпуска

Журнал вычислительной математики и математической физики

Архивы выпусков Информация о журнале Отправить рукопись в журнал