CC BY
348
УДК 519.245 Вестник СПбГУ. Сер. 1. 2013. Вып. 1
АЛГОРИТМ ПРИБЛИЖЕННОГО РЕШЕНИЯ ЗАДАЧИ КОММИВОЯЖЕРА
Т. М. Товстик1, Е. В. Жукова2
1. С.-Петербургский государственный университет, канд. физ.-мат. наук, доцент, peter.tovstik@mail.ru
2. С.-Петербургский государственный университет, студент, catich06@mail.ru
Введение. Коммивояжер должен найти минимальный замкнутый маршрут обхода N пунктов с посещением каждого из них ровно один раз. Это задача дискретной оптимизации NP-сложности, требующая для своего решения экспоненциально большого времени, так как нужно рассмотреть N!/(2N) вариантов. Иногда при больших N достаточно найти субоптимальное (близкое к оптимальному) решение.
За последние 30 лет размерности решаемых задач изменились от 100 до 108. Сначала, в основном, разрабатывались теоретические методы, а сейчас создаются также и комплексы программ, в которых решаются более общие задачи, частным случаем которых является задача о коммивояжере.
Метод под названием «сеть Хопфилда» был опубликован в 1985 г. [1]. Авторы статьи для решения задачи коммивояжера предложили использовать нейронные сети, минимизирующие энергию.
В 1987 г. Р. Дурбин и Д. Уиллшоу [2] предложили метод эластичной сети, который находил субоптимальное решение за меньшее время, чем с помощью сети Хопфилда. При N = 100 найденный вариант отличался от оптимального на 1%.
Метод имитации отжига, называемый чаще методом Метрополиса, был опубликован в 1953 г. [3]; он применяется для решения большого круга задач, частным случаем которых может служить задача о коммивояжере. На основе динамических методов Монте-Карло ищется глобальный минимум функционала.
Для поиска глобального экстремума с использованием метода Монте-Карло можно воспользоваться степенным методом, впервые упомянутым в монографии С. М.Ермакова [4] и в статье [5] описанным более подробно.
Появились модификации классической задачи коммивояжера. Сюда можно отнести задачу коммивояжера на максимум, заключающуюся в отыскании маршрута максимального веса с заданными весами дуг между пунктами. Асимптотически точный алгоритм можно найти в работе [6].
С приближенными алгоритмами решения задач на минимум и максимум для двух коммивояжеров можно ознакомится по работе [7].
В данной статье предлагается алгоритм построения хорошего начального приближения. Последующие приближения строятся с помощью моделирования динамических полей методом Метрополиса, который применяется на отдельных участках пути, что дает возможность воспользоваться параллельными вычислениями. Обсуждается выбор коэффициента отжига. Алгоритм приводит к маршруту, не имеющему самопересечений. Дается сравнение полученных результатов с результатами других
© Т. М. Товстик, Е. В. Жукова, 2013
авторов. Представленные примеры свидетельствуют о том, что предложенный алгоритм приводит к решениям, близким к оптимальным.
Известно [8], что при использовании метода Метрополиса к задачам с большим числом конфигураций важно начать процесс с хорошего начального приближения, так как именно оно будет уточняться при следующих итерациях. Рассмотрим вариант начального приближения, который затем будет уточняться методом Метрополиса.
Алгоритм построения начального приближения. Предположим, что заданы декартовы координаты х'г,у'г, г = , узлов (пунктов, точек). Путем линейных преобразований хг = х'о + Ххг, у[ = у'0 + Хуг переводим их в единичный квадрат. Введем полярную систему координат г,ф с началом в его центре. Точки нумеруем с ростом угла ф. Полученную нумерацию точек закрепляем и этот вариант считаем начальным приближением.
Последующие приближения основаны на использовании метода Метрополиса, в котором присутствует коэффициент отжига.
Рис. 1. Начальная нумерация точек (слева). Маршрут Hтz для 1000 случайных точек (справа).
Метод Метрополиса с коэффициентом отжига. Пусть X — конечное множество, его элементы х € X называются конфигурациями. Определяется вещественная функция Н(х), называемая энергетической, и равная длине пути в задаче коммивояжера. Метод Метрополиса дает возможность приближенно найти конфигурацию, минимизирующую энергетическую функцию Н(х).
В задаче коммивояжера заданы координаты N точек, конфигурация х определяется порядком их обхода. Так как маршрут замкнутый, номер, с которого начинается обход, не играет роли. Если на к-м приближении х = х(к), то
х« = -к - -к --- -к - = -к) , (^
.(к) .(к) где —номер пункта в начале маршрута, а -г —номер г-го пункта в порядке
обхода. При к = 0 ряд (1) определяет расположение точек в начальном приближении
х (0) = (1
—► 2 —► • •• —► N —► 1).
Энергетическая функция Н(х) принимается равной длине пути коммивояжера по замкнутому маршруту согласно конфигурации х = х(к):
N
Н (х) = Н (х(к)) = £ г(к), г(к) = г(3(к) ), (2)
г=1
I ■(к) -(к)Л
где т(з\ ,3т ) —расстояние между точками с соответствующими номерами.
Переход от к-го приближения к (к + 1)-му будем строить путем так называемой «двойной замены» (обращения), которая состоит в следующем. Случайно выбираются два номера пути обхода (1), пусть это оказались номера 3(к) и 3т), т > г. Формируется у — пробный вариант пути, в котором по сравнению с (1) частично меняется порядок обхода, а именно, города с номерами 3вк), г +1 ^ в ^ т — 1, обходятся в обратном порядке:
( (к) .(к) .(к) .(к) .(к) .(к) -(к) .(к)\ ,оч
у = (Я ) ^----► 3( ) ^ 3т-1 ^ 3т-2 ^ • • -3г(+1 ^ зт ^ 3т+1 ^----> 3^ ) • (3)
В алгоритме Метрополиса на (к + 1)-м шаге если Н(у) ^ Н(х(к)), то пробный вариант принимается в качестве нового приближения, т.е. х(к+1) = у. Если длина пробного пути Н(у) больше, чем у предыдущего приближения, т.е. Н(у) > Н(х(к)), то, чтобы не застрять в локальном минимуме, принимается х(к+1) = у с положительной вероятностью
Р = Р(х(к+1) = у)=ехр(-вДН), АН = Н (у) — Н (х(к)), (4)
а с вероятностью 1 — Р пробная конфигурация отклоняется и х(к+1) = х(к) •
Чтобы определить выбор конфигурации на (к + 1)-м шаге, с помощью датчика случайных чисел разыгрывается число а, равномерно распределенное на (0,1). Если оказывается, что а < Р, то х,(к+1) = у, иначе х,(к+1) = х(к •
В (4) величина в называется коэффициентом отжига [3].
Согласно теории моделирования динамических полей [8] при в ^ ж в пределе алгоритм Метрополиса приводит к конфигурации, минимизирующей длину пути коммивояжера.
О выборе коэффициента отжига в- Ясно, что коэффициент в должен быть соизмерим с ДН. Действительно, так как величины ДН могут сильно отличаться в разных задачах, при одинаковой вероятности принятия новой конфигурации выбор в будет зависеть от ДН. Выделим эту зависимость.
Из формул (1)—(4) находим
ДН = г(з(к) ,№_1) + ,з(к)) — г(з(к) ,3+1) — 3-1 ,з№). (5)
Сначала рассмотрим вариант со случайным расположением исходных точек. Будем предполагать, что точки равномерно распределены в единичном квадрате. В этом случае среднее значение и среднеквадратическое отклонение расстояния г между пунктами равны Ег = 0^5214, <г(г) = 0^2479 • Эти оценки получены методом Монте-Карло, и следовало бы использовать запись <г(г) « 0^2479, но здесь и в дальнейшем в таких случаях, подразумевая приближенный знак, будем писать знак равенства.
Поясним использование метода Монте-Карло в данном случае. Если координаты двух случайных точек равномерно распределены в единичном квадрате, то среднее значение Ег расстояния г между ними
Er
1111 V7(Ж1 - ж2)2 + (yi - y2)2dxidx2dyidy2-)
здесь все интегралы берутся в пределах от нуля до единицы. В методе Монте-Карло этот интеграл аппроксимируется суммой
1 п
Ег«гп = о = мл - «21?))2 + («3(я - «4(Я)2)1/2,
n j=i
в которой все ai(j), 1 ^ i ^ 4,1 ^ j ^ n, получаем с помощью датчика случайных чисел. С большой вероятностью можно утверждать, что \гп — Ег| = 0(а(г)/^/п). Здесь и в дальнейшем счет производился при n = 108. Подобным образом оценивались и другие величины, например, в формуле (11).
Представим в в виде
" = <6>
тогда при ДН > 0 вероятность (4) принятия пробной конфигурации y на (к + 1)-м шаге принимает вид
(7)
Методом Монте-Карло находим а(ДН|АН > 0) = 0.2946.
Обозначим через Po среднее значение вероятности (7) принятия новой конфигурации, когда АН > 0:
Связь во и Po представлена в таблице.
АН > 0 . (8)
/Зо 4 5 6 7 8 9 10 11 12 13
Ро 0.195 0.166 0.144 0.128 0.116 0.106 0.097 0.090 0.084 0.079
В дальнейшем метод Метрополиса будем применять к конфигурациям из ограниченных областей, и найденное ранее значение а(ДН) = 0.2945 не будет являться характеристикой ДН. Значит, нужно перейти к другому варианту задания в, который будет учитывать изменения длины маршрута при последующих итерациях.
Пусть для исходных данных разыгрываются точки для двойной замены, и при ДН > 0 отрезки £1 и £2 заменяются на П1 и П2. В результате получаем
ДН = П1 + П2 - £1 - £2 при ДН > 0. (9)
На к-й итерации формула (9) соответствует выражению (5). Дисперсию ДН представим в виде
а2(ДН|ДН > 0) = а2(т - £1) + а2(^ - £2) + а(т - £1 Мт - £2). (10)
Методом Монте-Карло оцениваем среднеквадратические отклонения а и коэффициенты корреляций р величин, входящих в состав ДН в формуле (9), при условии, что ДН > 0:
а(а) = а(Ь) = 0^2364, а(т) = а(п2) = 0-2328, а(т — £1) = а(п2 — £2) = 0-2883, Рт,^1 = Рт,Ь = °2449, Рт-£иП2-Ь = —04784 •
(11)
Из равенств (10) и (11) следует, что
а2(ДН\ДН > 0) « а2(щ — £1 |ДН > 0) = а2(т — Ь\ДН > 0) Аналогичные рассуждения приводят к результату
а2(П1 — £1) « 15а2 (£1) « 15а2(щ) « 075(а2(£1) + а2(т)). Отсюда получаем, что
а2(ДН\ДН > 0) « 0^75(а2(£1) + а2(щ))• Это дает возможность заменить формулу (6) на
(3= / 131 (12)
где /?! « fío/VOЛ5■
Полученные выше дисперсии и коэффициенты корреляций найдены для случайных независимых равномерно распределенных точек, а у маршрута после большого числа итераций будут другие параметры. Для маршрута после к-й итерации найдем среднее значение Г(к) расстояний между точками и разброс (а (к) )2 по формулам
^ = (*(к))2 = 4- (13)
Предлагается в качестве в использовать выражение вида (12),
/3= , (14)
в котором дисперсии заменены на (а(к))2 и (а(к+1))2, причем величина (а(к+1))2 вычислена у пробной конфигурации у.
Вероятность Р принятия в методе Метрополиса пробной конфигурации, если ДН > 0, будем вычислять по формуле
Р = ехр ( -/?« АИ | . (15)
Этой же формулой будем пользоваться, если конфигурация исходных точек задана и не является случайной.
Оптимизация маршрута на отдельных участках. При большом числе точек метод Метрополиса сходится медленно, поэтому для ускорения сходимости поступим
следующим образом. Будем стараться уменьшить общую длину за счет нахождения более коротких маршрутов на отдельных частях пути.
Назовем растром Г упорядоченную последовательность прохождения части пунктов
^ - - - (16)
размером растра — число пунктов Ь(Г) = т — г + 1 в ряде (16), а энергией растра — величину Н(Г), где
т— 1
н(Г) =53 гМ. (17)
в=1
Среди путей, соединяющих начальную и конечную точки растра и проходящих через все внутренние точки (16), методом Метрополиса ищем путь с энергией меньшей, чем Н(Г). При получении пробной конфигурации методом двойной замены первая и последняя точки в растре (16) не разыгрываются.
Указанным способом последовательно улучшаем части маршрута, пока не обойдем весь маршрут. У последовательно идущих растров должны быть непустые пересечения.
Сначала можно пройти по всему маршруту с размером растра порядка Ь « N/4, а затем, уменьшая или увеличивая размер растра, обойти маршрут несколько раз.
Если при осмотре полученного маршрута есть подозрение, что на каком-то участке маршрут не оптимальный, находим начало и конец растра, содержащего указанный участок, и применяем к нему метод Метрополиса.
Оценить полученный маршрут можно с помощью, так называемой, нормализованной длины маршрута
7 = Н/у/ЫА, (18)
где А — площадь области, содержащей все N точек. Известно [8], что при случайном независимом равномерном распределении точек внутри области выполняется неравенство 0.655 < ^ < 0.92 и, предположительно,
7 « 0.749. (19)
В принятом маршруте начального приближения нет самопересечений. Если на каком-то этапе оно возникло, его нужно убрать, причем общая длина пути при этом уменьшится. Действительно, путь АБМСП всегда длиннее пути АС МБП (рис. 2).
Рис. 2. Ликвидация самопересечений.
Расчеты показали, что самопересечения при длине пути Lo = L(BMC) < 10 можно, как правило, устранить, если пройти по всему маршруту с растром размера L, L0 + 1 < L < 15.
Итак, данный способ дает путь, близкий к оптимальному, а иногда и оптимальный. Этот вывод подтверждают приводимые ниже примеры.
При большом числе пунктов для ускорения решения задачи можно пользоваться параллельными вычислениями, оптимизирующими пути на отдельных участках.
Примеры. В первом примере исходные точки были промоделированы случайно, в остальных примерах данные взяты из библиотеки TSPLIB. У большинства примеров в библиотеке TSPLIB приводятся маршруты, которые были на момент их публикации оптимальными, и длины этих маршрутов будем обозначать Hopt. Длины маршрутов, полученные описанным выше способом, будем обозначать Htz.
Пример 1. Исходные данные являются реализацией 1000 независимых, равномерно распределенных в единичном квадрате точек. На левой части рис. 1 приведен путь обхода начального приближения, у которого H = 151.569,7 = 4.793. Справа на рис. 1 показан окончательный вариант маршрута с Htz = 23.929, 7 = 0.757, согласно приближенному равенству (19) он близок к оптимальному.
Рис. 3. Маршрут GR666.
Пример 2. Пример содержит координаты 666 пунктов. В библиотеке ТБРЫВ эти данные названы СИ666. В библиотеке под названием «оптимальный» приводится маршрут (изображен слева на рис. 3) с длиной Нр = 3952^53^ Он имеет несколько самопересечений, которые образуются от двух вертикальных линий. Длину пути можно уменьшить, если их убрать (вероятно, данные в библиотеке ТБРЫВ содержат опечатки). Маршрут, полученный описанным в статье способом, приведен справа на рис. 3. Он не имеет самопересечений, а его длина равна Hтz = 3240^46 и нормализованная длина 7 = 0^498, что меньше значения (19), ибо точки имеют сгущение. Область А в (18) — это прямоугольник, включающий все точки.
Пример 3. Пример СИ120 содержит данные 120 пунктов. На левом и правом рис.4 приведены пути, полученные соответственно М.Гретшелом в 1977 г. и описанным выше методом. В 2010 г. С.М.Ермаков и С.Н.Леора эти же данные обработали степенным методом, который описан в статье [5]. Там же приведен график полученного ими маршрута с длиной пути Неь = 1654^76, и она меньше, чем
Hgt = Hopt = 1666.51. С помощью алгоритма, приведенного в статье, удалось найти еще более короткий путь: Htz = 1645.60, 7 = 0.740. Все указанные маршруты не имеют самопересечений.
Рис. 4. Маршрут GR120.
Рис. 5. Маршрут PCB442.
Пример 4. Исходные данные РСВ442 взяты из библиотеки ТБРЫБ. Это задача Гретшела с 442 узлами, которая возникает при проектировании электронных плат. Приведенный в библиотеке оптимальный маршрут имеет длину Нр = 50.78. Для этого примера Hтz = 51.51, 7 = 0.716, т.е. полученный маршрут длиннее оптимального из библиотеки ТБРЫБ. На левом и правом рис.5 приводятся пути, соответствующие длинам Но^ и HTZ.
Пример 5. Пример содержит координаты 225 пунктов и в библиотеке ТБРЬ1Б имеет название Т8225. Оптимального варианта в библиотеке ТБРЬ1Б нет. На рис. 6 приводится наш вариант, не имеющий самопересечений. По всей видимости, это оптимальный вариант, хотя здесь возможны симметричные или круговые перестановки, не меняющие длину: Hтz = 126645.93, а 7 = 0.6978.
N Г 1
/ \ \ i А
V К ) \
1/1 N
Рис. 6. Оптимальный вариант для Т8225.
Приведем время (в секундах), которое компьютер затрачивает на получение приведенных в примерах маршрутов Hтz в зависимости от числа пунктов N:
Число пунктов 120 225 442 666 1000
Номер примера 3 5 4 2 1
Время в секундах 7 7 22 38 62
Таким образом, представленный алгоритм позволяет найти маршрут, близкий к оптимальному, за небольшое время работы компьютера. Особенностью алгоритма является то, что улучшение маршрута предполагает визуальное наблюдение за промежуточными результатами.
Литература
1. Hopfield J. J., Tank D. W. Neural computation of decisions in optimization problems // Biological Cybernetics. 1985. Vol.52. P. 141-152.
2. Durbin R., Willshaw D. An analogue approach to the traveling salesman problem using an elastic net method // Nature. 1987. Vol.326. 689 p.
3. Metropolis N., Rosenbluth A. W., Rosenbluth M. N., Teller A. H., Teller E. J. Equations of state calculations by fast computing machines //J. Chem. Phys. 1953. Vol.21. P. 1087-1092.
4. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. 472 с.
5. Ермаков С. М., Леора С. Н. Метод Метрополиса в задачах поиска экстремума // Математические модели. Теория и приложения. 2010. Вып. 11. С.-Петербург. C. 72-82.
6. Сердюков А. И. Задача коммивояжера на максимум в конечномерных вещественных пространствах // Дискретный анализ и исследование операций. 1995. Т. 27, №1. C. 50-56.
7. Гимади Э.Х., Глазков Ю.В., Глебов А.Н. Алгоритмы приближенного решения задачи о двух коммивояжерах в полном графе с весами ребер 1 и 2 // Дискретный анализ и исследование операций. 2007. Сер. 2. Т. 14, №2. C.41-61.
8. Винклер Г. Анализ изображений, случайные поля и динамические методы Монте-Карло. Новосибирск. Изд. СО РАН, 2002. 343 с.
9. Beardwood J., Haton J.H., Hammersley J.M. The shortest path through many points // Proc. Cambridge Phil. Soc. (1959). Vol.55. P. 299-327.
Статья поступила в редакцию 20 сентября 2012 г.
