В данной статье рассматриваем и реализуем на Python различные подходы (GA, ACO, SA) к решению задачи коммивояжёра.
В данной статье рассматриваем и реализуем на Python различные подходы (GA, ACO, SA) к решению задачи коммивояжёра.
О задаче
Задача коммивояжёра (Travelling Salesman Problem, TSP) — задача комбинаторной оптимизации. Как правило, её суть сводится к поиску оптимального пути, проходящего через все промежуточные пункты по одному разу и возвращающегося в исходную точку.
Для возможности применения математического аппарата к решению задачи мы представим её в виде симметричного графа. Вершины графа соответствуют условным городам, а рёбра — путям сообщения между этими городами. Каждому ребру сопоставим критерий выгодности маршрута — расстояние между городами. В целях упрощения задачи будем считать, что граф является полностью связным.
Для визуализации задачи будем использовать класс TSP (программный код представлен ниже). Этот класс включает в себя функции, необходимые для отрисовки точек, маршрутов, а также интерактивной легенды, на элементы которой можно нажимать для того, чтобы скрыть / отобразить выбранные пути. Сгенерируем набор точек, воспользовавшись функцией points = generate_problem(20), и выведем его на экран с помощью TSP(points=points, paths=None).
# Travelling Salesman Problem from random import randint from numpy import array from matplotlib import pyplot as plt from matplotlib.lines import Line2D from matplotlib.backend_bases import PickEvent from algorithms.utils.path import Path def generate_problem(count: int, canvas_size: int = 1000) -> list[tuple[int]]: """Generates a list of random 2D points.""" return [(randint(0, canvas_size), randint(0, canvas_size)) for _ in range(count)] class TSP: """ Allows to visualize the Traveling Salesman Problem and paths. """ CLR_POINT = "#eb343a" CLR_PATH = [ "#eb343a", "#db34eb", "#5b34eb", "#34b4eb", "#34eb4c", "#ebe534", "#eb9234", ] def __init__(self, points: list[tuple[int]], paths: list[Path] = None) -> None: """Initializes the problem, outputs its data using graphics.""" self.__points = points self.__paths = paths self.__fig, self.__ax = plt.subplots(num=f"Travelling Salesman Problem") self.__show() def get_points(self) -> list[tuple[int]]: """Getter to get the list of 2D points of the initialized problem.""" return self.__points def get_paths(self) -> list[Path]: """Getter to get the list of paths of the initialized problem.""" return self.__paths def __draw_points(self) -> None: """Draws 2D points and their coordinates on the canvas.""" self.__ax.scatter( *array(self.__points).T, zorder=1, color=TSP.CLR_POINT, label=f"Points ({len(self.__points)})", ) for i, p in enumerate(self.__points): plt.annotate( i + 1, p, ha="center", textcoords="offset points", xytext=(0, 4), fontsize=8, ) plt.annotate( f"({p[0]}; {p[1]})", p, ha="center", va="top", textcoords="offset points", xytext=(0, -4), fontsize=6, ) def __draw_paths(self) -> list[Line2D]: """Draws all given paths on the canvas.""" lines = [] if self.__paths: for i, path in enumerate(self.__paths): points = [self.__points[i] for i in path.indx] (l,) = plt.plot( *array(points).T, ls="--", zorder=0, color=TSP.CLR_PATH[i % len(TSP.CLR_PATH)], label=f"{path.name} ({path.leng:.2f})", ) lines.append(l) return lines def __draw_legend(self, lines: list[Line2D]) -> None: """Draws the legend on the canvas.""" if lines: self.__ax.set_title( "Tip: Click on the legend line(s) to turn the path ON / OFF", fontsize=10, loc="left", ) legend = self.__ax.legend() lined = {} for legline, origline in zip(legend.get_lines(), lines): legline.set_picker(True) lined[legline] = origline def on_pick(event: PickEvent) -> None: legline = event.artist origline = lined[legline] visible = not origline.get_visible() origline.set_visible(visible) legline.set_alpha(1.0 if visible else 0.2) self.__fig.canvas.draw() self.__fig.canvas.mpl_connect("pick_event", on_pick) else: self.__ax.legend() def __show(self) -> None: """Shows the canvas with the drawn data.""" self.__draw_points() lines = self.__draw_paths() self.__draw_legend(lines=lines) plt.show() if __name__ == "__main__": pass
Стоит отметить, что в программной реализации решениями является список объектов класса Path, которые содержат в себе информацию о последовательности точек маршрута, о его длине и названии (используются для легенды на графике).
# Path from dataclasses import dataclass @dataclass class Path: """ Dataclass describing a path using: * list of point indices; * path length; * path name (optional). """ indx: list[int] leng: float name: str
Но откуда берутся эти решения?
Поскольку коммивояжёр в каждом из условных городов встаёт перед выбором следующего города из тех, что он ещё не посетил, существует (n-1)! / 2 маршрутов для симметричной задачи. Если городов, например, 20, то возможных маршрутов будет 60 822 550 204 416 000. Для того, чтобы перебрать такое количество вариантов и найти оптимальный путь, понадобятся огромные вычислительные мощности и куча времени.
К счастью, на практике зачастую не требуется находить именно точное решение. Достаточно и субоптимального, которое будет близко к оптимальному пути.
При такой постановке задачи для её решения можно воспользоваться метаэвристическими алгоритмами, которые значительно ускорят процесс нахождения решения.
Метаэвристические алгоритмы
Метаэвриистика — метод оптимизации, многократно использующий простые правила для достижения субоптимального решения.
В данной статье мы рассмотрим следующие алгоритмы решения задачи:
Генетический алгоритм (Genetic Algorithm, GA), относящийся к классу эволюционных методов;
Муравьиный алгоритм (Ant Colony Optimization, ACO), относящийся к классу методов “роевого” интеллекта;
Алгоритм имитации отжига (Simulated annealing, SA), относящийся к классу методов, имитирующих физические процессы.
Каждый из вышеописанных алгоритмов будет использовать методы для вычислений расстояний между точками. Объединим данный функционал в класс Base, от которого в дальнейшем будем наследовать общие методы.
# Base from math import sqrt class Base: """ The base class for path finding algorithms. Contains common functions. """ @staticmethod def __euclidean_dist(a: tuple[int], b: tuple[int]) -> float: """Calculates the Euclidean distance between two 2D points.""" return sqrt((a[0] - b[0]) ** 2 + (a[1] - b[1]) ** 2) @staticmethod def _calculate_dist(dm: list[list[float]], indx: list[int]) -> float: """Calculates the path length based on the index list of the distance matrix.""" dist = 0 for i in range(len(indx) - 1): dist += dm[indx[i]][indx[i + 1]] return dist @staticmethod def _distance_matrix(points: list[tuple[int]]) -> list[list[float]]: """Calculates the distance matrix for the given 2D points.""" return [[Base.__euclidean_dist(a, b) for b in points] for a in points]
Здесь:
euclidean_dist() используется для нахождения Евклидова расстояния между двумя 2D точками;
calculate_dist() рассчитывает длину указанного пути на основе матрицы расстояний;
distance_matrix(), собственно, эту матрицу вычисляет.
Генетический алгоритм
Генетический алгоритм — алгоритм поиска, используемый для решения задач оптимизации и моделирования путём случайного подбора, комбинирования и вариации искомых параметров с использованием механизмов, аналогичных естественному отбору в природе.
Данные механизмы применяются к особям, каждая из которых в программной реализации представляет из себя список индексов 2D точек — последовательность, в которой эти точки следует обойти для достижения кратчайшего пути. Начальная популяция состоит из особей, случайно сгенерированных функцией initialization().
После создания начальной популяции запускается цикл операций, направленных на постепенную эволюцию особей:
fitness_sort() ранжирует особи в текущей популяции с помощью функции приспособления (чем короче маршрут, содержащийся в геноме особи, тем лучше);
selection() отбирает часть наиболее приспособленных особей из текущей популяции;
crossover() скрещивает особей текущей популяции, соединяя части генома от двух случайно выбранных родителей (панмиксия);
mutation() случайным образом делает перестановки в геноме некоторых особей текущей популяции.
Когда цикл завершит свою работу, мы используем лучшую особь (маршрут) из последней популяции в качестве решения задачи.
Отметим, что результаты работы алгоритма во многом зависят от заданных гиперпараметров:
population — общее число особей, участвующих в одной итерации;
iter — максимальное количество итераций алгоритма;
s — определяет, сколько лучших особей попадет в следующую популяцию (в процентном соотношении);
m — определяет, как часто особи в популяции мутируют (в процентном соотношении).
Ниже приводится программная реализация данного алгоритма.
# Genetic Algorithm from random import random, shuffle, sample from .utils.base import Base from .utils.path import Path class GA(Base): """ Genetic algorithm is a metaheuristic inspired by the process of natural selection that belongs to the larger class of evolutionary algorithms. ----- `population: int` THE NUMBER OF INDIVIDUALS The total number of individuals involved in one iteration. ----- `iter: int` THE NUMBER OF ITERATIONS The maximum number of iterations of the algorithm. ----- `s: float` SELECTION COEFFICIENT Determines how many of the best individuals will make it to the next population. ----- `m: float` MUTATION COEFFICIENT Determines how often individuals in a population mutate. """ def __init__(self, population: int, iter: int, s: float, m: float) -> None: """Initializes the hyperparameters for the algorithm.""" self.population = population self.iter = iter self.s = s self.m = m @staticmethod def __fitness_sort(dm: list[list[float]], individuals: list[list[int]]) -> None: """Sorts the individuals of a given population by fit.""" individuals.sort(key=lambda i: GA._calculate_dist(dm, i)) def __initialization(self, l: int) -> list[list[int]]: """Initializes the first population of individuals.""" base = list(range(l)) individuals = [] for _ in range(self.population): shuffle(base) individuals.append(base + [base[0]]) return individuals def __selection(self, individuals: list[list[int]]) -> None: """Selects the best individuals of a given population.""" del individuals[int(self.population * self.s) :] def __crossover(self, individuals: list[list[int]]) -> None: """Crossbreeding some individuals of a given population.""" childs = [] w_size = len(individuals[0]) // 2 for _ in range(len(individuals), self.population): p1, p2 = sample(individuals, 2) childs.append( p1[: w_size - 1] + [i for i in p2[:-1] if i not in p1[: w_size - 1]] + [p1[0]] ) individuals += childs def __mutation(self, individuals: list[list[int]]) -> None: """Mutates some individuals of a given population.""" sampling = list(range(1, len(individuals[0]) - 1)) for item in individuals: if random() < self.m: i, j = sample(sampling, 2) item[i], item[j] = item[j], item[i] def run(self, points: list[tuple[int]], name: str = None) -> Path: """Runs the algorithm for the given 2D points.""" l = len(points) dm = GA._distance_matrix(points) individuals = self.__initialization(l) for _ in range(self.iter): GA.__fitness_sort(dm, individuals) self.__selection(individuals) self.__crossover(individuals) self.__mutation(individuals) GA.__fitness_sort(dm, individuals) return Path(indx=individuals[0], leng=GA._calculate_dist(dm, individuals[0]), name=name) if __name__ == "__main__": pass
Такой результат получается при запуске алгоритма с набором гиперпараметров GA(population=1500, iter=40, s=0.2, m=0.5):
Муравьиный алгоритм
Муравьиный алгоритм — один из эффективных алгоритмов для нахождения приближённых решений задачи коммивояжёра, а также решения аналогичных задач поиска маршрутов на графах. Суть подхода заключается в анализе и использовании модели поведения муравьёв, ищущих пути от колонии к источнику питания.
Изначально условные муравьи передвигаются случайным образом, охватывая различные маршруты. Во время передвижения муравьи оставляют след из феромона, обозначая эти самые маршруты. Соответственно, чем короче получился путь, тем большее количество феромона за одинаковое количество времени муравьи на нём оставят. Феромон с течением времени испаряется. Из-за этого на длинных маршрутах остаётся меньше феромона, чем на коротких, даже если по ним передвигалось одинаковое количество муравьёв. Муравьи будут чаще выбирать самую “пахучую тропинку” и забывать пути с меньшим количеством феромона.
В программной реализации мы на каждой итерации цикла имитируем поведение каждого муравья из колонии:
create_indx() формирует маршрут муравья (список индексов точек) и подсчитывает длину;
select_i() отвечает за выбор муравьём каждой последующей точки маршрута внутри функции create_indx();
update_pm() обновляет матрицу феромона.
Формулой ниже описывается вероятность, с которой муравей перейдёт из точки i в точку j (используется в функции create_indx()):
Здесь:
?(i,j) — количество феромона на ребре i, j (значение из матрицы феромона);
? — гиперпараметр, контролирующий влияние ?(i,j) (задаётся пользователем);
?(i,j) — привлекательность ребра i, j, равняется 1 / d(i,j) (обратному значению из матрицы расстояний);
? — гиперпараметр, контролирующий влияние ?(i,j) (задаётся пользователем).
После формирования маршрутов необходимо обновить матрицу феромона (за это отвечает функция update_pm()). Обновление происходит по следующему правилу:
Здесь:
?(i,j) — количество феромона на ребре i, j (значение из матрицы феромона);
? — скорость испарения феромона (задаётся пользователем);
??(i,j) — количество отложенного феромона, равняется q / d(i,j) (интенсивность феромона, делённая на длину маршрута из точки i в точку j).
В конце каждой итерации находим лучший маршрут из маршрутов, пройденных на текущей итерации, и сравниваем его с абсолютно лучшим найденным маршрутом. Если новое решение оказалось лучше, то записываем его в отдельную переменную. Последний записанный маршрут и будет являться решением задачи.
Так же, как и генетический алгоритм, данный алгоритм зависим от гиперпараметров:
ants — общее количество агентов (муравьев), задействованных в одной итерации;
iter — максимальное количество итераций алгоритма;
? — коэффициент, контролирующий влияние феромона на ребре;
? —коэффициент, контролирующий влияние привлекательности маршрута (величина, обратно пропорциональная длине маршрута);
? — коэффициент испарения феромона, который отражает степень взаимного влияния между муравьями, как правило, значение равно [0, 1], что предотвращает бесконечное накопление феромона;
q — интенсивность феромона, которая представляет собой общее количество феромонов, в определённой степени влияет на скорость сходимости алгоритма.
Ниже приводится программная реализация данного алгоритма.
# Ant Colony Optimization from math import inf from random import random, shuffle from .utils.base import Base from .utils.path import Path class ACO(Base): """ Ant Colony Optimization algorithm is introduced based on the foraging behavior of an ant for seeking a path between their colony and source food. ----- `ants: int` THE NUMBER OF ANTS The total number of agents (ants) involved in one iteration. ----- `iter: int` THE NUMBER OF ITERATIONS The maximum number of iterations of the algorithm. ----- `a: float` INFORMATION ELICITATION FACTOR The information elicitation factor ?, which represents the relative importance of the pheromone, reflects the importance of the accumulation of the pheromone with regard to the ants' path selection. ----- `b: float` EXPECTED HEURISTIC FACTOR The expected heuristic factor ?, which represents the relative importance of the visibility, reflects the importance of the heuristic information with regard to the ants' path selection. ----- `p: float` PHEROMONE EVAPORATION COEFFICIENT The pheromone evaporation coefficient ?, which represents the degree of pheromone evaporation, reflects the degree of mutual influence among ants. Generally, the value of is [0, 1], which prevents the infinite accumulation of pheromone effectively. ----- `q: float` PHEROMONE INTENSITY The pheromone intensity Q, which represents the total pheromone, affects the convergence speed of the alghoritm to a certain extent. """ def __init__(self, ants: int, iter: int, a: float, b: float, p: float, q: float) -> None: """Initializes the hyperparameters for the algorithm.""" self.ants = ants self.iter = iter self.a = a self.b = b self.p = p self.q = q @staticmethod def __select_i(selection: list[int]) -> int: """Selects a random index of the next 2D point.""" sum_num = sum(selection) if sum_num == 0: return len(selection) - 1 tmp_num = random() prob = 0 for i in range(len(selection)): prob += selection[i] / sum_num if prob >= tmp_num: return i def __create_indx(self, dm: list[list[float]], pm: list[list[float]]) -> list[int]: """Creates a new ordering of 2D point indices based on the distance and pheromone.""" l = len(dm) unvisited_indx = list(range(l)) shuffle(unvisited_indx) visited_indx = [unvisited_indx.pop()] for _ in range(l - 1): i = visited_indx[-1] selection = [] for j in unvisited_indx: selection.append( (pm[i][j] ** self.a) * ((1 / max(dm[i][j], 10**-5)) ** self.b) ) selected_i = ACO.__select_i(selection) visited_indx.append(unvisited_indx.pop(selected_i)) visited_indx.append(visited_indx[0]) return visited_indx def update_pm(self, pm: list[list[float]], tmp_indx: list[list[int]], tmp_leng: list[float]) -> None: """Updates the pheromone matrix.""" l = len(pm) for i in range(l): for j in range(i, l): pm[i][j] *= 1 - self.p pm[j][i] *= 1 - self.p for i in range(self.ants): delta = self.q / tmp_leng[i] indx = tmp_indx[i] for j in range(l): pm[indx[j]][indx[j + 1]] += delta pm[indx[j + 1]][indx[j]] += delta def run(self, points: list[tuple[int]], name: str = None) -> Path: """Runs the algorithm for the given 2D points.""" l = len(points) dm = ACO._distance_matrix(points) pm = [[1 for _ in range(l)] for _ in range(l)] res_indx = [] res_leng = inf for _ in range(self.iter): tmp_indx = [] tmp_leng = [] for _ in range(self.ants): indx = self.__create_indx(dm, pm) tmp_indx.append(indx) tmp_leng.append(ACO._calculate_dist(dm, indx)) self.update_pm(pm, tmp_indx, tmp_leng) best_leng = min(tmp_leng) if best_leng < res_leng: res_leng = best_leng res_indx = tmp_indx[tmp_leng.index(best_leng)] return Path(indx=res_indx, leng=res_leng, name=name) if __name__ == "__main__": pass
Такой результат получается при запуске алгоритма с набором гиперпараметров ACO(ants=100, iter=20, a=1.5, b=1.2, p=0.6, q=10):
Алгоритм имитации отжига
Алгоритм имитации отжига — метод решения задачи глобальной оптимизации. Он основывается на имитации физического процесса, который происходит при кристаллизации вещества, в том числе при отжиге металлов. Предполагается, что в определённый момент атомы вещества практически выстроены в кристаллическую решётку, но ещё допустимы переходы отдельных атомов из одной ячейки в другую. Активность атомов тем больше, чем выше температура, которую постепенно понижают, что приводит к тому, что вероятность переходов в состояния с большей энергией уменьшается. Устойчивая кристаллическая решётка соответствует минимуму энергии атомов, поэтому атом либо переходит в состояние с меньшим уровнем энергии, либо остаётся на месте.
В программной реализации под энергией понимается длина получившегося маршрута. Изначально маршрут задаётся случайно. На каждой итерации с помощью перестановки пары индексов мы генерируем новый маршрут. На основе его длины, длины предыдущего маршрута и текущего значения температуры функция is_acceptable() определяет, совершать ли переход в новое состояние или нет. Вероятность перехода рассчитывается по формуле:
Здесь:
?F(i) — изменение “энергии” (разность между длинами вероятного и текущего маршрута);
T(i) — значение температуры на i-ой итерации алгоритма.
В конце каждой итерации система “остывает”: параметр, отвечающий за температуру, уменьшается. Когда алгоритм отработает отведённое количество итераций, он вернёт минимальный найденный маршрут, являющийся решением задачи.
Гиперпараметры алгоритма:
iter — максимальное количество итераций алгоритма;
t — начальная температура поиска, уменьшается по мере продвижения поиска;
g — коэффициент, влияющий на изменение температуры.
Ниже приводится программная реализация данного алгоритма.
# Simulated Annealing # from math import exp from numpy import exp from random import sample, random from .utils.base import Base from .utils.path import Path class SA(Base): """ Simulated annealing is a probabilistic technique for approximating the global optimum of a given function. Specifically, it is a metaheuristic to approximate global optimization in a large search space for an optimization problem. ----- `iter: int` THE NUMBER OF ITERATIONS The maximum number of iterations of the algorithm. ----- `t: int` INITIAL TEMPERATURE The initial temperature for the search decreases with the progress of the search. ----- `g: float` CHANGE COEFFICIENT The coefficient affecting temperature change. """ def __init__(self, iter: int, t: int, g: float) -> None: """Initializes the hyperparameters for the algorithm.""" self.iter = iter self.t = t self.g = g def __is_acceptable(self, prb_leng: float, tmp_leng: float) -> bool: """Checks if the state transition will execute.""" prob = min(1, exp(-(prb_leng - tmp_leng) / self.t)) if prob > random(): return True return False def run(self, points: list[tuple[int]], name: str = None) -> Path: """Runs the algorithm for the given 2D points.""" l = len(points) dm = SA._distance_matrix(points) tmp_indx = [i for i in range(l)] + [0] tmp_leng = SA._calculate_dist(dm, tmp_indx) res_indx = tmp_indx.copy() res_leng = tmp_leng for _ in range(self.iter): i, j = sample(range(1, l), 2) prb_indx = tmp_indx.copy() prb_indx[i], prb_indx[j] = prb_indx[j], prb_indx[i] prb_leng = SA._calculate_dist(dm, prb_indx) if self.__is_acceptable(prb_leng, tmp_leng): tmp_indx = prb_indx tmp_leng = prb_leng if tmp_leng < res_leng: res_indx = tmp_indx res_leng = tmp_leng self.t *= self.g return Path(indx=res_indx, leng=res_leng, name=name) if __name__ == "__main__": pass
Такой результат получается при запуске алгоритма с набором гиперпараметров SA(iter=20000, t=100, g=0.6):
Сравнение алгоритмов
Сравнивать алгоритмы затруднительно, потому что — как упоминалось выше — результаты во многом зависят от подобранных гиперпараметров. Так же нужно учитывать, что существует довольно много различных вариаций реализации алгоритмов, в зависимости от которых результаты тоже будут отличаться.
Если же выбирать лучшую реализацию из тех трёх, что были описаны в данной статье, то я бы отдал предпочтение муравьиному алгоритму. Он довольно быстр и, что наиболее важно, показывает хорошие результаты даже при большом количестве исходных данных. Ниже приводятся примеры решения задачи на 100 и 200 точек.
Заключение
Таким образом, мы разобрали одни из самых популярных метаэвристических подходов к решению задачи коммивояжёра: генетический алгоритм, муравьиный алгоритм и алгоритм иммитации отжига.
Подробнее ознакомиться с кодом данных алгоритмов и поэкспериментировать с подбором гиперпараметров можно на GitHub’е по данной ссылке.