CC BY
18
УДК 519.863
М. О. Солнцева1, Б. Г. Кухаренко1'2
1 Московский физико-технический институт (государственный университет) Институт машиноведения РАН
Применение методов кластеризации узлов на графах с разреженными матрицами смежности в задачах
логистики
Рассматривается решение задачи оптимального расположения логистических центров в сети поставок на графе транспортной сети с помощью различных алгоритмов кластеризации. Разбиение узлов транспортной сети на кластеры выполняется после предварительного анализа данных, представляемых сильно разреженной матрицей смежности графа этой сети. В качестве альтернативного подхода к кластеризации узлов на графе рассматриваются процедуры построения и усечения минимального дерева Штейнера.
Ключевые слова: анализ данных, анализ графов, разреженные матрицы смежности, алгоритмы кластеризации, деревья Штейнера, транспортные сети, логистика.
Введение
При проектировании логистических сетей решается задача оптимального расположения логистических центров [1]. Принято выделять три основные модели расположения логистических центров: создание единого центра по обслуживанию потребителей (single sink network design problem); расположение нескольких новых центров (multi sink network design problem); разбиение всех потребительских центров на подмножества и создание в каждом подмножестве своей сети (selective network design problem) [2]. Точное решение этих задач на неориентированном графе принадлежит к классу NP-сложных задач.
Настоящая работа посвящена решению задачи оптимального расположения логистических центров при проектировании цепи поставок (facility location problem) с помощью алгоритма кластеризации K-means [3] и его модификаций [4] и также с помощью построения и усечения дерева Штейнера [5]. При решении этой задачи существует ограничение, налагаемое топологией транспортной сети, которая задаётся матрицей смежности (affinitv matrix) графа этой сети. Для построения графов транспортной сети используются данные проекта OpenStreetMap [6]. На рис. 1 приведён пример графа транспортной сети Омского региона. До настоящего времени ограничения, налагаемые реальной транспортной сетью, при решении задачи оптимизации не принимались во внимание. Введение такого ограничения необходимо для регионов с транспортной сетью, пред ставимой разреженной матрицей смежности. При решении поставленной задачи разбиение графа транспортной сети на произвольное число кластеров нецелесообразно. Необходимо предварительно выявить скрытую структуру графа, которая отражена в его матрице смежности. Анализ матрицы смежности позволяет оценить количество кластеров в структуре графа. Такую оценку можно использовать в качестве входного параметра для алгоритмов кластеризации узлов, заданных на графе.
Алгоритм ALA (Alternate Location Allocation) [7], обычно используемый для решения этой задачи, не учитывает топологию транспортной сети, поэтому его результаты применимы только для областей с хорошо развитой транспортной инфраструктурой [8]. Существенным недостатком этого алгоритма является то, что его целевая функция не достигает глобального минимума и зависит от начальных условий, задаваемых случайным образом. Наиболее близким к алгоритму ALA является алгоритм кластеризации K-means [3].
Рис. 1. Граф транспортной сети Омского региона, построенный по данным OperiStreetMap
Для того чтобы оценить максимально допустимое число логистических центров, в настоящей работе проводится предварительный анализ скрытой структуры графа транспортной сети, представленного разреженной матрицей смежности. Такой анализ выполняется с помощью алгоритма кросс-ассоциаций [9], выделяющих) в матрице смежности кластеры с однородной внутренней структурой, в результате чего определяется число кластеров в структуре графа, которое является входным параметром для применяемых далее алгоритмов кластеризации. Алгоритмы K-means [3] и K-means++ [4| используются для решения задачи выделения нескольких логистических центров в сети поставок (multi sink network design problem). Процедуры построения и усечения дерева Штейнера [5] применяются для выделения как нескольких логистических центров (multi sink network design problem), так и единого центра (single sink network design problem).
В настоящей работе разработан комплекс программ на языке С++ и других объектно-ориентированных языках, реализующих указанные выше алгоритмы кластеризации [3, 4, 5, 9], который позволяет производить численные эксперименты на графах реальных транспортных сетей.
Материалы и методы
Данные OSM
Для построения графа транспортной сети используются массивы данных проекта OpenStreetMap [6]. В настоящей работе исследуется граф транспортной сети Омского региона (см. рис. 1). Исходный граф включает 140 ООО узлов. Это число сокращается до 38 777 за счет исключения несущественной информации. В описании графа сохраняются узлы, соответствующие исходному пункту каждой дороги, пункту ее окончания и пересечениям с другими дорогами. Объем данных, представляющий эту разреженную матрицу смежности с 100 378 ненулевых элементов, занимает на диске 3 Гб и 1,5Гб в памяти компьютера.
OMSK REGION'S ROAD MAP
57.5
71.5 72.0 72.5 73.0 73.5 74.0 74.5 75.0
EL,0
Предварительный анализ матрицы смежности графа транспортной сети. Алгоритм кросс-ассоциаций
В настоящей работе для предварительного анализа, выявляющего скрытую структуру графа, используется метод кросс-ассоциаций [9]. Этот метод предлагает общую модель, основанную на сжатии данных согласно принципу MDL (Minimum description length) [10].
Структура графа отображается в его матрице смежности - бинарной матрице из нулей и единиц, поэтому проблема её анализа аналогична задачам анализа ассоциативных правил, информационного поиска и выделения групп в сетях и т.п. [9]. Во всех этих случаях маркировка столбцов и строк не имеет значения. Группировка независимых строк и столбцов матрицы основывается на их сходстве. Эти перекрёстные ассоциации выявляют скрытую структуру графа. Соответствующие прямоугольные области различной плотности используются для быстрого продвижения по структуре матрицы.
Если исходные элементы матрицы составляют m, х п «прямоугольников» с «плотностью» либо 0, либо 1, то по завершении процедуры выявляются прямоугольные блоки с плотностью от 0 до 1, число которых подлежит определению. Отбрасывание части этих прямоугольников уменьшает сложность описания данных. Для оценки сложности описания матрицы смежности в используемом алгоритме кросс-ассоциаций [9] применяется MDL-принцип [10], в котором стоимость структуры данных, выделяемой в матрице, оценивается числом бит, необходимых для передачи всей структуры вместе с данными о каждой выделенной прямоугольной области.
Пусть D = [dij ] — матрица бинарных данных размерности т, х п, где тшп количество строк и столбцов матрицы соответственно. Припишем строки к группам строк и столбцы к группам столбцов:
Ф : {1, 2,...,т} ^ {1, 2,...,к}, Ф : {1, 2,...,п} ^ {1, 2,..,1},
где к и I обозначают число независимых групп среди строк и столбцов соответственно. Согласно этому приписанию (или «перекрёстной ассоциации» {Ф, Ф}) элементы матрицы D перестраиваются таким образом, чтобы строки, соответствующие группам 1, 2 и т.д., перечислялись в соответствии с этим порядком и аналогично - для столбцов. Такая перестройка разделяет исходную матрицу D на прямоугольные блоки меньшего размера -подматрицы Dij размерное ти (a,i,bj ), г = 1, 2,... ,к, j = 1, 2,... ,1.
Введём обозначения: А - бинарная матрица, размерности а х b, п i(A) - число ненулевых элементов матрицы А, п о(А) - число нулевых элементов матрицы А, п(А) = п i(A) + по(А) = а х b, РА(г) = ^ , i = 1, 0. Тогда полная длина кода в битах
С (А) = £ пг(А) log ( = п(А)Н (РА(0)),
где Н - бинарная энтропия Шэннона. Полная длина кода для матрицы D с учётом заданной кросс-ассоциации имеет вид
Т(D; к;, I; Ф; Ф) =
1-1 1-1 к I к I
= log* к + log* I + ^ [logOi\ + ^ [logb3] + ^ ^ [logM, + 1)1 + ^ ^С(А,,), 3 = 1 j = 1 i=1 3 = 1 i=1 3 = 1
где Гц = (Y<k=i ^ - к + г, г = 1,..., к - 1-, Ь, = (Т,[=з k) - I + 3, 3 = 1,..., I - 1.
Оптимальная кросс-ассоциация соответствует числу групп строк к*, числу групп столбцов I* и кросс-ассоциации {Ф*, Ф*} таких, что полная результирующая длина кода Т (D; к*; 1*;Ф*;Ф*) минимизируется.
Задача определения оптимальной кросс-ассоциации является вычислительно сложной, поэтому для её решения применяется эвристика [9], состоящая из двух шагов:
1) Внутренний цикл: для заданных чисел к и I найти подходящую перегруппировку (т.е. кросс-ассоциацию), соответствующую достижению локального минимума функции
к I
ЕЕ с )• а)
г=1 з=1
2) Внешний цикл: поиск наилучших к и I среди тех, которые рассматриваются во внутреннем цикле. Этот этап использует резкое падение функции (1) при малых значениях к и I (см. рис. 2, где к = I = И).
Рис. 2. Значения функции (1) в зависимости от выделенного числа однородных кластеров в матрице смежности графа транспортной сети Омской области (см. рис. 3)
о
1-4
G
N
Рис. 3. Результат применения алгоритма кросс-ассоциаций к матрице смежности графа транспортной сети Омской области, представленного на рис. 1
Оценка сложности алгоритма: на каждом шаге алгоритма возрастает либо либо /, поэтому сумма к + I всегда возрастает на 1. Следовательно, общая сложность алгоритма
оценивается как 0(щ(0)(к* + I*)2. На практике для нахождения кросс-ассоциации достаточно 20 итераций.
На рис. 3 показаны результаты применения алгоритма поиска кросс-ассоциаций [9] для матрицы смежности графа транспортной сети Омской области.
Алгоритм кластеризации К-теат
Пусть (XI,..., Хм} — векторное представление кластеризуемых узлов, {ц^} — набор векторов, соответствующих центрам кластеров, гпк € (0,1} - индикаторные переменные, указывающие на приписание узлов кластерам. Алгоритм К-теапв [3] минимизирует суммарное среднеквадратичное отклонение узлов кластера от его центра:
т1п {У ] У Упк ||х„ - ^к У2
{N К }
^^rnk ||Xn-ßkII2 ,
п=1к=1 )
откуда ßk = ^-k'ГпкХ" .
У rt ^nk
Сложность алгоритма определяется как 0(N 2К+1 log N), N - число узлов, К - число кластеров. Этот алгоритм гарантирует нахождение локальной оптимальной конфигурации кластеров, но не применим для нахождения глобального оптимума. Алгоритм K-means [4] имеет по меньшей мере два основных недостатка. Первый из них - время работы алго-
К
_ Пр0изв0льная аппроксимация недостаточно хорошо оптимизирует целевую функцию. Чтобы избежать главным образом второго недостатка, используется алгоритм K-means++ [4]. Сначала этот алгоритм осуществляет процедуру оптимальной инициализации центров кластеров, затем выполняются стандартные итерации алгоритма K-means [3]. Благодаря процедуре инициализации сложность алгоритма K-means++ определяется как 0(log К).
Алгоритм кластеризации K-means++
Рассмотрим n-мерное евклидово пространство (хотя подход может быть представлен и в более общем виде). Пусть имеется множество точек X и требуется найти множество С, состоящее из К точек, называемых центрами кластеризации, которые минимизируют ошибку
Erг(С, X) = -L. ^d(x,C(х)),
1 1 х£Х
С( x) С x
С (х) = argmind(x, с).
На первом шаге алгоритма происходит выбор первого центра с [1] кластеризации из однородного распределения точек множества X. Далее выполняется выбор k-ro центра с[к] ( k = 2,..., К) из мультиномиального распределения на множестве X, где точка х имеет вероятность
*х = Ёхда - «<х)2.
D(x) определяется как расстояние до ближайшего существующего центра.
После того как начальные центры выбраны, переходят к стандартной процедуре кластеризации K-means [3].
Процедуры, построения и усечения дерева Штейнера
Алгоритм построения и усечения дерева Штейнера также можно использовать для выделения центров кластеризации. Полное описание алгоритмов построения дерева Штейнера (Steiner tree) можно найти в работе [5].
Пусть G = (V, Е) — неориентированный граф, в котором каждое ребро ij £ Е имеет неотрицательную реальную стоимость, и пусть N — множество концевых узлов (вершин) графа, N С V. Дерево Т, выделяемое в графе G, называют деревом Штейнера, если оно содержит все концевые узлы графа V(Т) 5 N. Требуется найти такое дерево Штейнера Т, стоимость которого с(Т) = ^ijeE(T) C(Ü) минимальна среди всех деревьев Штейнера для графа G. Оптимальное решение Т является минимальным деревом Штейнера.
Для усечения дерева Штейнера строят сеть расстояний Dq(N), где N — выделенное указанное множество концевых узлов, и затем находят минимальное остовное дерево М в Dq(N). Заменяя каждое ребро v w в М кратчайшим путём v — w в G, получают новый подграф Т ' в G. В результате находят минимальное остовное дерево Т', последовательно удаляя любые узлы Штейнера первой степени.
Результаты и обсуждения
Рис. 4. Схематичное представление структуры матрицы смежности графа транспортной сети Омской области (см. рис. 3). Чем темнее прямоугольник, тем выше плотность связей между узлами в нём
В настоящей работе задача оптимизации расположения логистических центров в условиях ограничений решается с помощью алгоритмов кластеризации после того, как выявлена структура разреженной матрицы смежности графа рассматриваемой транспортной сети. Следует отметить, что ранее проблема оптимизации расположения логистических
центров при ограничениях транспортной сетью в литературе не рассматривалась.
Результаты предварительного анализа графа транспортной сети Омской области (см. рис. 3), выполненные с помощью метода кросс-ассоциаций [9], показали, что в структуре этого графа выделяются 15 кластеров, расположенных на главной диагонали матрицы смежности графа. На рис. 4 показано идеализированное изображение структуры матрицы смежности, представленной на рис. 3. Насыщенность цвета прямоугольников на рис. 4 зависит от величины плотности связей в каждом выделенном кластере. Чем темнее прямоугольник, тем больше плотность связей в кластере.
Отметим, что полученное на этом этапе число кластеров отражает всего лишь общую структуру матрицы смежности, в то время как для размещения логистических центров имеет смысл рассматривать только те кластеры, плотность связей в которых выше некоторого среднего уровня. На рис. 5 показано распределение значений плотности связей в выделенных 15 кластерах. Как видно из рис. 5, только 8 кластеров имеют плотность связей выше среднего значения. Это число кластеров следует рассматривать как предельное количество логистических центров, к которым можно приписать выделенные узлы графа. В зависимости от числа рассматриваемых узлов можно варьировать число размещаемых логистических центров.
х 10""
Cluster densities
0.9
0.8
0.7
O.B
0.5
0.4
0.3
0.2
0.1
# cluster densities mean density
* I
i 1
1 1 i 1
3:56C Я'00'6 ► 1 »
î-04 ►
Ф
i ► < ►
^
ft 4 ► 1 ►
1 2 3 4 5 В 7
9 10 11 12 13 14 15
cluster index
Рис. 5. Точки графика показывают значения плотности связей в каждом из 15 выделенных кластеров. Сплошная линия соответствует средней плотности связей в выделенных кластерах
Для решения задачи оптимального расположения новых логистических центров в сети поставок в настоящей работе рассматривались наиболее популярные алгоритмы класте-
ризации, такие как К-теапв и его модификации [3, 4], и также процедуры построения и усечения дерева Штейнера [5]. Исходные целевые функции этих алгоритмов не учитывают естественные ограничения, возникающие в конкретных приложениях, где наиболее значительным ограничением является транспортная сеть. В настоящей работе такие ограничения были учтены. Для построения графа транспортной сети были использованы массивы данных проекта Ореп81:тее1,Мар [6]. Работа программ демонстрируется на примере региона Омской области.
На рис. 6 приведены примеры кластеризации с помощью алгоритмов К-теапв [3] и К-теапвН—Ь [4]. При этом искомые центры кластеризации соединяются с каждым рассматриваемым узлом по маршруту, вычисляемому алгоритмом Дейкстры [11]. С помощью алгоритма К-теапв-Ы- [4] можно не только существенно уменьшить время работы алгоритма К-теапв [3], но и получить более рациональное разбиение на кластеры, соответствующее минимальному расстоянию между узлами в кластере, как это видно из рис. 6 (б).
Рис. 6. Сравнение результатов кластеризации алгоритма К-теапв (а) и К-теапв-Н- (б) в условиях ограничений, налагаемых транспортной сетью, на карте Омской области. Центры кластеров обозначены кружками большего диаметра
Рис. 7. Результаты работы алгоритма построения минимального дерева Штейнера (а) и процедуры усечения этого дерева для выделения трёх кластеров (б)
Построение дерева Штейнера на неориентированном графе с заданными п концевыми узлами начинается с построения сети расстояний (distance network) [5]. Положение кон-
цевых узлов задаётся пользователем на карте дорог. В соответствии с алгоритмом среди множества деревьев Штейнера для заданных концевых узлов находится минимальное дерево. Результаты алгоритма, реализующего построение минимального дерева Штейнера, представлены на рис. 7а.
Полученное дерево представляет собой минимальный маршрут, соединяющий все заданные узлы. Любой узел может быть выделен в качестве логистического центра, если решается задача выбора расположения единого центра в сети (single sink network design). Если ставится задача выбора нескольких логистических центров, то к полученному минимальному дереву Штейнера (рис. 7а) сначала применяется процедура его усечения. На рис. 76 показано усечение минимального дерева Штейнера с целью выделения трёх кластеров. Усечение производится последовательным удалением наиболее длинных рёбер полученного дерева.
Заключение
Предварительный анализ структуры транспортной сети, выполненный на матрице смежности её графа, показал, что число кластеров, выделенное в структуре матрицы смежности, может быть использовано как оценка предельного значения числа логистических центров.
Литература
1. Facility Location: Applications And Theory / Ed. by Drezner Z., Hamacher H.W. — BerlinHeidelberg: Springer-Verlag, 2004.
2. Awerbuch В., Azar Y. Buv-at-bulk network design //In IEEE Symposium on Foundations of Computer Science (FOCS). - 1997. - P. 542-547.
3. MacQueen J. Some Methods for Classification and Analysis of Multivariate Observations // Proc. Fifth Berkeley Svmp. Math. Statistics and Probability. — 1967. — P. 281-296.
4. Arthur D., Vassilvitskii S. K-means++: the advantages of careful seeding / SODA'07 Proceedings of the Eighteenth Annual ACM-SIAM Symposium on Discrete Algorithms.
— CitvPhiladelphia, StatePA: placecountrv-regionSIAM Press. — 2007. — P. 1027-1035.
5. Cheriyan J., Ravi R. Lecture Notes on Approximation Algorithms for Network Problems
— Canada: University of Waterloo, 1998.
6. Ramm F., Topf J., Chilton S. OpenStreetMap: Using and Enhancing the Free Map of the World // Cambridge, United Kingdom: UIT Cambridge Ltd. — 2010. — P. 386.
7. Cooper L. Location-allocation problems // Operations Research. — 1963. — V. 11. — P. 331 343.
8. Michael G.K., Ajithkumar N.P. Material Flow Analysis of Public Logistics Networks // Progress in Material Handling Research. — 2002. — P. 205-218.
9. Chakrabarti D., Papadimitriou S., Modha D., Faloutsos C. Fully automatic cross-associations // Proceedings of the Tenth ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. — 2004. — P. 79-88.
10. Rissanen J. Information and Complexity in Statistical Modeling // Springer. — 2007. — P. 97-103.
11. Dijkstra E.W. A note on two problems in connection with graphs // Numerische Mathematik. - 1959. - I. 1. - P. 269-271.
Поступила в редакцию 14-01.2012.
