CC BY
23
Головешкин В.А.1, Жукова Г.Н.2, Ульянов М.В.3, Фомичев М.И.4
1 Институт прикладной механики РАН, Московский государственный университет информационных технологий, радиотехники и электроники (МИРЭА), г. Москва, д.т.н., профессор,
nikshevolog@yandex.ru
2 Московский государственный университет печати (МГУП) им. Ивана Федорова, г. Москва,
к. ф.-м. н., доцент, gnzh@mail.ru
3 Институт проблем управления РАН им. В.А. Трапезникова, г. Москва, д.т.н., профессор,
muljanov@mail.ru
4 ФКН НИУ Высшая школа экономики (ФКН НИУ ВШЭ), г. Москва, студент, michan94@yandex.ru
СРАВНЕНИЕ РЕСУРСНЫХ ХАРАКТЕРИСТИК ТРАДИЦИОННОГО И МОДИФИЦИРОВАННОГО МЕТОДА ВЕТВЕЙ И ГРАНИЦ ДЛЯ TSP7
КЛЮЧЕВЫЕ СЛОВА
Задача коммивояжера (TSP), метод ветвей и границ (МВГ), ресурсная эффективность, время работы алгоритма, бета распределение II типа, система функций Пирсона.
АННОТАЦИЯ
Сравниваются ресурсные характеристики модифицированного и классического МВГ для TSP. На основании экспериментальных результатов показано, что по величине затраченного на поиск решения времени модифицированный вариант МВГ эффективнее классического. Исследована стохастическая зависимость между временем работы каждого из двух исследуемых вариантов МВГ при фиксированном порядке матрицы стоимостей. Также описана зависимость затраченной памяти и времени работы алгоритма от порядка матрицы стоимости.
Постановка задачи
Задача коммивояжера относится к ЫР-трудным задачам комбинаторной оптимизации и заключается в отыскании оптимального (например, самого дешевого) гамильтонова цикла в полном ориентированном графе. Граф задается матрицей стоимостей, внедиагональные элементы которой (неотрицательные числа) интерпретируются как расстояния между пунктами назначения или стоимости перевозок. Диагональные элементы матрицы стоимостей полагаются условно равными бесконечности, чтобы осуществить запрет петель.
Впервые метод ветвей и границ (МВГ) для решения задачи коммивояжера был предложен в [1], хотя предварительные результаты и общие формулировки идей были представлены в более ранних работах [2-6]. Хорошо исследована задача коммивояжёра с симметричной матрицей, изучен асимметричный случай при условии, что выполняется неравенство треугольника (т.е. перевозка напрямую из А в Б выгоднее, чем через один или несколько промежуточных пунктов).
В дальнейшем для решения задачи коммивояжера применялись эволюционные, в основном генетические и муравьиные алгоритмы, а также комбинации метода ветвей и границ с генетическими алгоритмами [7]. Генетические и муравьиные алгоритмы использовались для повышения точности оценки снизу, что позволило более эффективно исключать бесперспективные подзадачи.
Несмотря на то, что задача коммивояжера достаточно хорошо исследована, прогнозирование трудоемкости или временных характеристик для конкретной задачи вызывает существенные трудности. Проводились исследования по прогнозированию времени, требующегося для решения задачи коммивояжера методом ветвей и границ, на основе случайного выбора [9]. В [10] предложена процедура прогнозирования числа вершин дерева решений на основе анализа поддерева решений, полученного на начальном этапе расчетов. Особенностью перечисленных методов прогноза является зависимость прогноза от конкретного экземпляра задачи, а также то, что прогноз может быть получен только спустя какое-то время после начала расчетов.
7 Набота выполнена при поддержке РФФИ, грант №16-07-160
Особенности модифицированной реализации МВГ
Поиск решения задачи коммивояжера методом ветвей и границ (МВГ) заключается в разделении множества допустимых решений на подмножества с целью дальнейшего сокращения перебора. На каждом шаге выбирается дуга ветвления, и множество туров делится на классы, один из которых состоит из туров, содержащих эту дугу, а второй — из оставшихся допустимых туров (ветвление). В каждом классе производится оценка минимальной стоимости тура и находится стоимость какого-нибудь допустимого тура из этого класса. Если стоимость найденного тура меньше стоимости ранее найденного «рекордного» тура, то новый тур объявляется «рекордом». Классы допустимых туров с нижней оценкой, превышающей стоимость «рекорда», исключаются из дальнейшего рассмотрения как бесперспективные. Если стоимость «рекорда» такова, что все нижние оценки других классов больше, то найден оптимальный тур. Каждый класс допустимых туров подвергается процедуре ветвления, пока не будет найдено оптимальное решение.
Метод ветвей и границ состоит в построении и исследовании древовидной структуры пространства решений задачи коммивояжера. Построение поискового дерева решений начинается с корня, который будет соответствовать множеству «всех возможных туров», т. е. корень дерева представляет множество R всех (п — 1)! возможных туров в задаче с п городами. Ветви, выходящие из корня, отличаются наличием или отсутствием в турах одной дуги, например, (к, /). Множество R делится на два класса { к,! } (туры из R , содержащие дугу (к,!) ) и { к,! } (туры, не содержащие эту дугу). Подробное поэтапное изложение метода можно найти, например, в [2,3].
Пример фрагмента такого дерева приведен на рисунке 1 (источник рисунка [2]).
Рис. 1. Фрагмент поискового дерева решений
Вычислительный эксперимент
Генерация матриц стоимости несимметричной задачи коммивояжера проводилась стандартным для языка С++ псевдослучайным равномерным генератором. Элементы матриц (размерности от 11 до 46 с шагом 1) были целыми числами из диапазона от 1 до 215 . Был реализован алгоритм [1] (без предвычисления начального тура жадным алгоритмом) со структурой хранения поискового дерева решений в виде бинарной кучи в классической реализации с перевычислением матрицы стоимостей (т.е. без использования дополнительной памяти) и в модифицированной реализации с хранением локальных (для вершин поискового дерева) матриц стоимостей (далее будем для краткости называть эту модифицированную реализацию МВГ1).
Для каждого фиксированного n выполнялось 10 000 псевдослучайных генераций матрицы стоимостей порядка n, каждый элемент матрицы - псевдослучайное число с дискретным равномерным распределением в диапазоне от 1 до 215 ■ Для каждой матрицы была решена задача коммивояжера классическим МВГ (с перевычислением матриц) и МВГ1 (с хранением матриц). Для каждого запуска программы замерялось время, затраченное на поиск оптимального тура, зависящее от числа вершин порожденного поискового дерева решений, а для МВГ1 также и максимальное значение затраченной на хранение матриц дополнительной памяти.
Эксперименты проводились на стационарном компьютере со следующими характеристиками:
— процессор: Intel i7 3770K 3800 ГГц
— оперативная память: Kingston KHX1600C9D3P1 16 ГБ
— материнская плата: GIGABYTE GA-Z77X-D3H
— операционная система: Fedora 21 workstation
Алгоритмы реализованы на языке С++ (компилятор gcc 4.9.2 20150212, Red Hat 4.9.2-6). Взаимозависимость временных характеристик модифицированного и классического
МВГ
Обозначим t1 время, затраченное классическим МВГ на поиск оптимального тура коммивояжера, t2 - время работы МВГ1, т2 - объем памяти, использованный для хранения промежуточных матриц. Изобразим результаты вычислительного эксперимента на диаграмме рассеяния (рис. 2)
Рис. 2. Диаграмма рассеяния t1 и t2 , n = 20, 30,40
На рис. 2 изображена синяя линия регрессии, подобранная по методу наименьших квадратов. Отметим, что с ростом размерности матрицы стоимостей уменьшается коэффициент наклона этой прямой, что свидетельствует об увеличении «превосходства» модифицированного метода над классическим в аспекте экономии затраченного времени. Зависимость коэффициента наклона slope линии регрессии t1 на t2 от размерности n матрицы стоимостей (а также коэффициента детерминации r2 ) изображена на рис. 3.
Рис. 3. Зависимость slope иг от n
Близкие к единице значения коэффициента детерминации указывают на тесную линейную зависимость t1 и t2 , что свидетельствует о большей экономичности (в смысле затрат времени) МВГ1.
Характеристики распределения времени работы МВГ1
Время t , затраченное на поиск оптимального тура коммивояжера, будем рассматривать как случайную величину, определенную на множестве всевозможных матриц размерности n с элементами — натуральными числами, не превосходящими 215 . Это пространство конечное, но содержит огромное число элементов, поэтому для определенной выше случайной величины будем подбирать такое непрерывное распределение, что его функция распределения не слишком сильно отличается от выборочной функции распределения полученных данных.
Данные t1 представляют собой положительные числа, равные времени, затраченному алгоритмом МВГ на решение задачи коммивояжера размерности n , t2 — время работы МВГ1
над теми же задачами, т2 — объем памяти, затраченной МВГ1 на хранение матриц. Задача решалась для п от 11 до 46, в каждом случае было проведено 10 000 вычислительных экспериментов.
Сначала проанализируем гистограммы результатов вычислительного эксперимента. На рис. 4, 5 и 6 приведены гистограммы времени t1 работы МВГ, t2 (времени работы МВГ1) и т2 — объема памяти, затраченной МВГ1 на хранение матриц. На каждом рисунке три гистограммы, соответствующие матрицам стоимости порядка п= 20, 30 , 40 .
Рис. 4. Гистограммы t1 , п= 20, 30 , 40
Рис. 5. Гистограммы t2 , п= 20, 30 , 40
Рис. 6. Гистограммы т2 , п= 20, 30 , 40
Все гистограммы похожи по форме, что и следовало ожидать, учитывая тесную линейную зависимость t1 и t2 (аналогичная зависимость имеет место и для т2 и t2). Вид гистограмм побуждает подбирать возможное распределение экспериментальных данных в классе непрерывных распределений с плотностью, представляющей собой убывающую на [а;Ь ] функцию, или же имеющую один максимум вблизи а , где а>0 Ь<ю .
Будем искать подходящее распределение, имеющее конечный четвертый момент. Для таких распределений определены коэффициент асимметрии Y1 и коэффициент эксцесса Y2 . Пусть наблюдается случайная величина X , имеющая математическое ожидание ЕХ , дисперсию DX , третий центральный момент у3 (X ) = Е (X — ЕХ )3 и четвертый центральный момент у4 (X ) = Е (X — ЕХ )4 , тогда коэффициент асимметрии Y1 и коэффициент эксцесса Y2 определяются по формулам [11]
У1 =
X)
Е ( X—EX)
(DX )3/2 (Е (X—EX )2)3
= у 4 (X) = Е (X—EX )2
Y 2=( DX )2 =( Е (X—EX )2 )2.
3
Покажем, что у 1 и у2 инвариантны относительно линейных преобразований X , для этого преобразованием сдвига и масштаба построим по X случайную величину Y :
Y=kX+x0 . (1)
Выразим моменты и коэффициенты асимметрии и эксцесса Y через соответствующие характеристики X .
1.Матожидание EY=E (kX+x0) =kEX+x0 .
2.Дисперсия
DY=E(Y-EY )2 =Е (kX+x0-(kEX+x0) )2 =Е(kX-kEX )2 =к2 DX .
3.Третий центральный момент E (Y- EY )3 =E ( kX- kEX )3 =k3E (X-EX )3 .
4.Четвертый центральный момент E (Y-EY )4 =k4 E (X-EX )4 .
5.Коэффициент асимметрии
(Y)_ E(Y- EY )3 = k3E(X- EX )3 = k3E (X-EX )3 = E (X-EX )3
™ Ь (DY )3/2 "(k2 E (X—EX )2 Г" k3 (E (X—EX )2)"2" (E (X - EX )2)3/211 .
6.Коэффициент эксцесса
E (Y - EY )4 _ k4 E (X - EX )4 _ E (X - EX )4 _
У 2 (Y ) = ^--3 =-±-^-3= /-V -3 = У2 ( X )=У 2 .
(DY) ^ (DX) (DX)
Как видно, коэффициенты асимметрии и эксцесса не изменяются при линейном преобразовании, что позволяет подбирать вид распределения по выборочным коэффициентам асимметрии и эксцесса, а затем за счет подбора коэффициентов k и X0 добиваться равенства выборочных среднего и дисперсиии Y теоретическим значениям для подобранного распределения.
В соответствии с подобранным распределением вычисляются математическое ожидание EY и дисперсия DY , затем по ним находятся k и X0 из (1):
k=J— , хо =EY — kEX . (2)
\DX
Рис.7. Коэффициенты асимметрии и эксцесса ( Г2 )
Следует отметить, что по значениям выборочных коэффициентов асимметрии и эксцесса не всегда возможно подобрать хотя бы одно «классическое» распределение.
Для подбора вида распределения воспользуемся методом К. Пирсона [13]. На рис. 7 в системе координат «коэффициент асимметрии - коэффициент эксцесса» изображены некоторые
классические распределения, а также точки, соответствующие значениям коэффициентов асимметрии и эксцесса экспериментальных данных. Данные представляют собой 10000 положительных чисел, равных времени t2 , затраченному модифицированным алгоритмом МВГ на решение задачи коммивояжера размерности п . Задача решалась для п от 11 до 46, каждому из этих случаев соответствует фиолетовый маркер.
Выборочные коэффициенты асимметрии и эксцесса t1 (время работы классического МВГ), t2 (время работы модифицированного МВГ), а также т2 (объем затраченной на хранение матриц памяти) изображаются на рисунке 8 точками примерно из одной и той же области.
к-1 аснмыегрим
Рис.8. Коэффициенты асимметрии и эксцесса ( t1 , t2 , m2 )
Если результаты эксперимента возможно описать единообразно некоторым распределением, то это распределение должно иметь коэффициенты асимметрии и эксцесса «между» соответствующими значениями для гамма и логнормального распределения, это область распределений шестого типа системы Пирсона (бета распределение второго типа, Beta prime distribution) [12,13].
Бета распределение второго типа описывается двумя параметрами а> 0 и в> 0 , причем четвертый момент конечен только при в> 4 . Плотность распределения имеет вид [12]:
а—11 1 X ( 1 + X
где
f (Х )=: В (а, в) • В (а, в) — бета функция, коэффициенты асимметрии и эксцесса равны соответственно
2(2а+в-1) I в-2
г 1
У 2 = 6-
(1 в—3 » а (а+в—1) = 6а (а+в—1)(5в—11) + (в—1 )2(в—2) = а (а+в—1)(в—3)(в—4)
= 6(5 в—11) +6( в—1 )2 (в—2) =B+C_
(в—3)(в—4) а(а+в—1)(в—3)(в—4) а(а+в—1)'
(3)
(4)
(5)
В_6(5в-11) г _6 (в-1 )2(в-2 ) В_(в-3)(в-4) 'С_(в-3)(в-4) •
Параметр ( находится из (5) как решение квадратного уравнения относительно (
а(а+в-1 )= С „=1-^, (6)
У2- О 2
где D = ( в-1) +-„ . Перед дискриминантом выбран знак плюс, чтобы обеспечить
У2-°
положительное значение а с учетом в> 4 .
и D
в-1
Подставляя С и D в (6) получим
а =
24(в-2 )_+ 1-1
(7)
ГУ2 (в-3)(в-4 )-6 (5 в-11) /
Затем возведем обе части (4) в квадрат (обе части положительны, поскольку в> 4 )
2= 4 (2 а + в-1 )2 (в-2) У1 а (а+в-1 )(в-3 )2
и придадим полученному уравнению вид
а(а+в-1)(в-3)2у2-4(2а+в-1 )2(в-2)=0 . (8)
Параметр в найдем, подставив в (8) а из (7) и решив полученное уравнение:
12+2 2 в = 3 +-.
2У2-3У1
Подставляя найденное значение в в (7), найдем а .
Некоторые значения параметров, полученные по экспериментальным данным приведены в таблице 1 и на рис. 9.
Табл. 1. Точечные оценки параметров бета распределения 11 типа
5 12 т2
п а в а в а в
21 0.656233 14.78672 0.771754 23.20291 0.851123 19.771012890
22 0.259704 5.144625 0.672774 5.228451 1.196278 5.443517978
23 0.626576 9.684045 0.72864 15.68017 0.842763 13.592175809
24 0.373771 6.602246 0.606926 6.985486 0.68882 7.012765161
25 0.427045 8.047167 0.594325 9.015214 0.620072 8.944139370
26 0.528095 9.268592 0.819729 9.190158 0.871106 9.230012948
27 0.416076 15.87569 0.692569 9.277104 0.755955 8.554410508
28 0.493696 14.66053 0.659561 12.06239 0.642138 13.137767265
29 0.273568 6.170698 0.270182 5.814129 0.443837 6.279142455
30 0.370631 5.464604 0.420186 5.485285 0.524322 6.653685070
31 0.25838 10.17488 0.256808 8.034325 0.353448 7.460506640
32 0.380466 7.469743 0.339419 7.011105 0.459047 8.019133906
33 0.271364 8.240574 0.254145 8.05216 0.350575 8.009574496
34 0.117049 7.045287 0.100416 7.475167 0.165922 6.909838078
35 0.236292 5.401035 0.212038 6.056703 0.299458 6.838542935
36 0.052658 5.186972 0.060503 5.204344 0.086996 5.065607549
37 0.069879 7.094084 0.06939 7.53807 0.113203 7.665997073
38 0.253724 13.82903 0.213711 16.2068 0.272672 12.530710979
39 0.195116 14.98218 0.171445 15.1198 0.213455 11.395406361
40 0.102776 13.48028 0.111546 12.99175 0.154693 9.911065410
41 0.112186 6.277624 0.126788 6.524766 0.174313 7.811211803
42 0.110775 7.706199 0.123706 7.295158 0.160005 6.678007491
43 0.123975 16.2609 0.130429 18.62152 0.161407 12.659731130
44 0.115105 10.82112 0.125164 9.71375 0.153199 8.375463681
45 0.124168 6.660459 0.145357 8.200421 0.165672 10.185727312
46 0.124124 7.273358 0.150179 7.998998 0.188199 11.502324733
MIN 0.052658 5.144625 0.060503 5.204344 0.086996 5.065607549
MAX 0.656233 16.2609 0.819729 23.20291 1.196278 19.77101289
(dpl.il
Рис.9. Точечные оценки параметров бета распределения 11 типа
На рис. 10 приведены гистограмма для t1 при п_ 40 и плотность бета распределения II типа с параметрами, рассчитанными по выборке.
Рис.10. Гистограмма t1 и плотность бета распределения II типа
Выводы
Модифицированный вариант МВГ по показателю «затраченное на поиск решения время» эффективнее классического варианта, причем выигрыш во времени (для п_ 40 почти в пять раз) становится все более значимым с ростом размерности матрицы стоимостей. Между значениями затраченного времени обоих вариантов МВГ наблюдается тесная линейная зависимость при фиксированном п , теснота связи растет с ростом п , приближаясь к 1.
Выигрыш во времени достигается ценой значительных затрат памяти, которые также растут с ростом п . Данные вычислительных экспериментов достаточно хорошо аппроксимируются бета распределением II типа, выбранным по предложенной методике на основе коэффициентов асимметрии и эксцесса.
Литература
1. Little, J. D. C., K. G. Murty, D.W. Sweeney, C. Karel. 1963. An algorithm for the traveling salesman problem. Operations Research 11, 972-989.
2. Dantzig, G., R. Fulkerson, S. Johnson. 1954. Solution of a large scale traveling salesman problem. Technical Report P-510. RAND Corporation, Santa Monica, California, USA.
3. Dantzig, G. B., D. R. Fulkerson, S. M. Johnson. 1959. On a linearprogramming, combinatorial approach to the traveling-salesman problem. Operations Research 7, 58-66.
4. Eastman, W. L. 1958. Linear Programming with Pattern Constraints. Ph.D. Thesis. Department of Economics, Harvard University, Cambridge, Massachusetts, USA.
5. Land, A. H., Doig A. G. 1960. An automatic method of solving discrete programming problems. Econometrica 28, 497-520.
6. Manne, A. S., H. M. Markowitz. 1956. On the solution of discrete programming problems. Technical Report P-711. RAND Corporation, Santa Monica, California, USA.
7. C. Cotta, J. Aldana, A. Nebro, J. Troya, Hybridizing genetic algorithms with branch and bound techniques for the resolution of the TSP, in: D. Pearson, N. Steele, R. Albrecht (Eds.), Artificial Neural Nets and Genetic Algorithms 2, Springer-Verlag, Wien New York, 1995, pp. 277-280.
8. G. Carpaneto, P. Toth, Some new branching and bounding criteria for the asymmetric traveling salesman problem, Management Science. 1980, v. 26, 736-743.
9. Knuth, D. E. 1975. Estimating the efficiency of backtracking programs. Mathematics of Computing, v. 29, pp. 121-136.
10. Cornuejols, G.; Karamanov, M.; and Li, Y. 2006. Early estimates of the size of branch-and-bound trees. INFORMS Journal on Computing. 2006, v. 18, No. 1, pp. 86—96.
11. Крамер Г. Математические методы статистики. М.: Мир, 1975. — 648 с.
12. Jonhnson, N.L., Kotz, S., Balakrishnan, N. 1995. Continuous Univariate Distributions, v. 2, Wiley.
13. Pearson, K Contributions to the Mathematical Theory of Evolution. III. Regression, Heredity and Panmixia, Philosophical Transactions of the Royal Society of London, 1896. 187, pp. 253-318.
14. Ульянов М.В., Фомичев М.И. Ресурсные характеристики способов организации дерева решений в методе ветвей и границ для задачи коммивояжера // Бизнес - информатика, 2015, (статья принята к публикации в №4 журнала в 2015 г.)
15. Ульянов М.В. Ресурсно-эффективные компьютерные алгоритмы. Разработка и анализ. — М.: ФИЗМАТЛИТ, 2008. — 304 с.
