Все пути одинаковы: они ведут в никуда. Но у одних есть сердце, а у других — нет. Один путь дает тебе силы, другой — уничтожает тебя.
- Карлос Кастанеда
Здравствуйте уважаемые дамы и господа, а также не бинарные личности. Хорошей эпохи.
Мы с вами уже пробовали решать точно задачу коммивояжёра методом динамического программирования и методом ветвей и границ, результат не плох, но слабоват.
В данной статье постараюсь показать, что точное решение ближе, чем принято считать.
Мы будем использовать метод целочисленного программирования, который является частным случаем линейного программирования, который в свою очередь является подклассом математического программирования.
Математическое программирование - это область математики, разрабатывающая теорию, численные методы решения многомерных задач оптимизации с ограничениями.
Ограничения могут быть:
дискретными (дискретное программирование или комбинаторная оптимизация (метод динамического программирования, ветвей и границ можно отнести сюда))
линейными (целочисленными или непрерывными)
нелинейными
Но и сама целевая функция может быть:
линейной (Linear programming, LP)
квадратичной (Quadratic programming, QP)
другой нелинейной (Conic optimization, др.)
В данной работе нас будет интересовать линейное программирование, где целевая функция и ограничения выражены линейно.
В свою очередь от того какой тип переменных используется, можно выделить:
Булево программирование (Binary integer programming, BIP)
Целочисленное (все переменные целочисленные, Integer programming, IP)
Смешанное (часть переменных целочисленная, Mixed-integer programming, MIP)
Линейное (все переменные рассматриваются как не целые, непрерывные – continuous, LP)
Для каждого класса задач есть свои методы, часть пересекается, часть имеет свои границы применения. Я не буду заострять ваше внимание на этом. Меня в первую очередь интересует прикладная применимость в работе.
Поэтому, если все сильно упростить, то можно сделать вывод: чем более непрерывны переменные, тем легче найти решение (Linear programming легче находит решение, чем Mixed-integer и Integer programming)
Но так как в нашей жизни многие процессы дискретны (например, трудно представить покупку 1,25 автомобиля, или разбить 0,5 яйца), на практике применяется метод релаксации (аппроксимация сложной проблемы соседней проблемой, которую легче решить).
То есть, задача решается как не целочисленная (или не бинарная), а затем к полученным данным применяются дополнительные ограничения, которые приводят задачу к желаемому результату.
Прежде чем мы перейдём непосредственно к алгоритму, скажем пару слов за инструментарий.
Решать задачу я предлагаю на языке Python, как очень удобном средстве прототипирования и проектирования с большим набором внешних библиотек расширения.
Для наших целей кроме основного набора модулей нам понадобится решатель (solver) задач целочисленного (IP) или смешанного программирования (MIP). Подойдёт любой корректный решатель.
В своих изысканиях я опробовал три библиотеки содержащих решатели IP/MIP: PuLP, CVXPY и SciPy.
PuLP - даже из названия становится понятно, что это библиотека, которая предназначена для решения задач линейного программирования. Среди иных альтернатив она выделяется минимализмом, проста для освоения новичками, отличается удобством использования и возможностью подключать внешние решатели.
CVXPY - это встроенный в Python язык моделирования с открытым исходным кодом для задач выпуклой оптимизации. В отличие от PuLP, может решать не только линейные задачи. Содержит более богатый функционал и приличный набор подключаемых решателей. Однако, как мне показалось CVXPY справляется с однотипными задачами медленнее, чем PuLP на тех же внешних решателях (по крайней мере, MIP) и более сложен для изучения. Если вам не требуется для работы выпуклая целевая функция, то я бы его для данной задачи не рекомендовал.
SciPy - библиотека для языка программирования Python с открытым исходным кодом, предназначенная для выполнения научных и инженерных расчётов. Нас интересует её модуль optimize – средства оптимизации. Скажу сразу, это весьма сложный подход, но среди всех опенсорсных решений самый быстрый, что я нашёл.
В примерах к статье буду приводить примеры в основном с использованием PuLP, чтобы вам легче было разобраться.
Ниже будут рассмотрены существующие решатели только целочисленных и смешанных линейных уравнений, применимых для наших целей.
Сначала установим сам PuLP он собирается со встроенным решателем CBC от COIN-OR Foundation, имеет неплохую производительность.
pip install pulp
SCIP - так же свободный решатель, но для данной задачи он меня несколько разочаровал. Это не означает что решатель плох, как по мне он уступает CBC в скорости.
CPLEX – классный коммерческий решатель от IBM. Существует комьюнити версия не для коммерческого использования, которая может решать линейные уравнения до 1000 переменных и 1000 ограничений (для задачи коммивояжера матрицы n<=45)
pip install cplex docplex
GUROBI - ещё один замечательный коммерческий решатель. Существует комьюнити версия не для коммерческого использования, которая может решать линейные уравнения до 2000 переменных (для задачи коммивояжера матрицы n<=63)
pip install gurobipy
Xpress - так же хороший проприетарный решатель.
Комьюнити версия не для коммерческого использования, может решать линейные уравнения до 5000 переменных + ограничений (для задачи коммивояжера матрицы ~n<=98)
pip install xpress
Здесь не просто так приведены ссылки на коммерческие решатели, так как они справляются с той же задачей в несколько раз быстрее на одинаковом наборе данных. У CPLEX и GUROBI есть возможность скачать бесплатную академическую лицензию. Так как я не являюсь участником университетской деятельности, мне этот путь заказан, но согласитесь, подход компаний подкупает, пока не посмотришь на коммерческие ценники.
Перейдём к самому алгоритму.
Шаг
Возьмём симметричную матрицу смежности для некого графа на 5 вершин:
[[ 0 14 3 13 8]
[14 0 16 2 7]
[ 3 16 0 14 10]
[13 2 14 0 6]
[ 8 7 10 6 0]]
Найдём для неё минимальный Гамильтонов цикл используя PuLP.
import pulp as pl
# Создадим модель задачи для решения - минимизировать длину пути
model = pl.LpProblem(name="tsp", sense=pl.LpMinimize)
# Подключаем решатель
solver = pl.PULP_CBC_CMD(msg=False)
# Заводим массив переменных
x = [pl.LpVariable(name=f'x_{i:02}_{j:02}', cat='Binary') for i in range(n) for j in range(n)]
# Устанавливаем целевую функцию
model += pl.lpSum([matrix[i][j] * x[i * n + j if i < j else j * n + i] for i in range(n) for j in range(n) if i < j])
# Вносим ограничения модели
for i in range(n):
model += pl.lpSum([x[i * n + j if i < j else j * n + i] for j in range(n) if i != j]) == 2
Если сейчас вывести модель, то она будет выглядеть так:
tsp:
MINIMIZE
14*x_00_01 + 3*x_00_02 + 13*x_00_03 + 8*x_00_04 + 16*x_01_02 + 2*x_01_03 + 7*x_01_04 + 14*x_02_03 + 10*x_02_04 + 6*x_03_04 + 0
SUBJECT TO
_C1: x_00_01 + x_00_02 + x_00_03 + x_00_04 = 2
_C2: x_00_01 + x_01_02 + x_01_03 + x_01_04 = 2
_C3: x_00_02 + x_01_02 + x_02_03 + x_02_04 = 2
_C4: x_00_03 + x_01_03 + x_02_03 + x_03_04 = 2
_C5: x_00_04 + x_01_04 + x_02_04 + x_03_04 = 2
VARIABLES
0 <= x_00_01 <= 1 Integer
0 <= x_00_02 <= 1 Integer
0 <= x_00_03 <= 1 Integer
0 <= x_00_04 <= 1 Integer
0 <= x_01_02 <= 1 Integer
0 <= x_01_03 <= 1 Integer
0 <= x_01_04 <= 1 Integer
0 <= x_02_03 <= 1 Integer
0 <= x_02_04 <= 1 Integer
0 <= x_03_04 <= 1 Integer
x_[i]_[j] – это переменная которая отвечает за то, выбирается данное ребро или нет (1 / 0).
14 * x_[i]_[j] – это вес ребра из матрицы смежности.
Ограничения говорят нам о том, что в каждую вершину должно вести ровно два пути, по одному мы приходим в вершину, а по-другому можем выйти.
Запустив решатель
status = model.solve(solver)
for v in model.variables():
print(f'{v.name} = {round(v.value())}')
print(f'Minimum path = {round(model.objective.value())}')
получим:
x_00_01 = 0
x_00_02 = 1
x_00_03 = 0
x_00_04 = 1
x_01_02 = 0
x_01_03 = 1
x_01_04 = 1
x_02_03 = 1
x_02_04 = 0
x_03_04 = 0
Minimum path = 34
Нам остаётся только отобразить на графе данные рёбра.
Мы нашли Гамильтонов цикл минимальной длины.
Шаг
Возьмём матрицу крупнее и проведём те же действия
[[ 0 5 8 5 7 1]
[ 5 0 12 4 4 5]
[ 8 12 0 13 15 7]
[ 5 4 13 0 2 6]
[ 7 4 15 2 0 8]
[ 1 5 7 6 8 0]]
О нет, что-то пошло не так!
Minimum path = 26
Мы нашли минимальный путь, проходящий через все вершины, но вот граф распался на два множества, так дело не пойдёт.
Как мы понимаем, что решение не верное? У нас есть верный признак: если множество распадается на два или более подмножества, значит это неверное решение.
В модуле networkx есть масса всяческих полезных алгоритмов для работы с графами на все случаи жизни.
Нас будет интересовать функция nx.connected_components, которая возвращает генератор для графа в виде набора множеств, из которого он состоит.
print(list(nx.connected_components(graph))
[{0, 2, 5}, {1, 3, 4}]
Нам нужно ввести дискриминирующую функцию, которая бы отвергала решение в случае, если оно распадается на множества. Мы поступим следующим образом: если в графе больше одного множества, то мы откажемся от такого решения и добавим новое ограничение. Ограничение должно быть таким, что бы мешало сформировать ошибочное решение повторно.
model += pl.lpSum([x[i[0] * n + i[1] if i[0] < i[1] else i[1] * n + i[0]] for i in graph.edges]) <= n - 2
В модели сформировалось новое условие:
_C7: x_00_02 + x_00_05 + x_01_03 + x_01_04 + x_02_05 + x_03_04 <= 4
Запускаем решатель ещё раз, с уточнёнными ограничениями.
Minimum path = 31
Путь стал немного длиннее чем ранее, но зато решение верно.
Шаг
Для графов большей размерности мы можем проворачивать подобный трюк в цикле до тех пор, пока решение не будет найдено.
Однако даже для такого не серьёзного графа, как этот, n=13:
[[ 0 35 87 67 77 63 51 118 121 132 116 9 92]
[ 35 0 52 58 45 30 21 91 109 118 107 44 73]
[ 87 52 0 71 20 39 39 57 98 103 104 97 61]
[ 67 58 71 0 51 80 44 63 54 65 50 74 30]
[ 77 45 20 51 0 45 26 48 81 87 85 86 42]
[ 63 30 39 80 45 0 36 92 123 131 125 71 85]
[ 51 21 39 44 26 36 0 70 90 99 90 60 53]
[118 91 57 63 48 92 70 0 54 53 64 127 33]
[121 109 98 54 81 123 90 54 0 12 13 128 38]
[132 118 103 65 87 131 99 53 12 0 24 140 46]
[116 107 104 50 85 125 90 64 13 24 0 122 43]
[ 9 44 97 74 86 71 60 127 128 140 122 0 100]
[ 92 73 61 30 42 85 53 33 38 46 43 100 0]]
Алгоритму потребовалось пройти 208 шагов и добавить 207 дополнительных ограничений, к тому же расчёт затянулся на несколько минут.
Так как с ростом размера входной матрицы количество множеств, на которое может распасться решение возрастает экспоненциально. Нам нужна другая стратегия!
Шаг
Утверждаю, что аналогичного решения можно добиться всего за 4 шага решателя и добавив всего 6 ограничений. Всё что нужно было сделать, это изменить набор ограничивающих условий.
Дискриминирующая функция должна создать такие условия, что бы препятствовали распаду графа на отдельные множества.
Вот одно из ограничений:
Зелёным обозначены рёбра которые попадают в неравенство
_C14: x_00_02 + x_00_03 + x_00_04 + x_00_05 + x_00_06 + x_00_07 + x_00_08
+ x_00_09 + x_00_10 + x_00_12 + x_01_02 + x_01_03 + x_01_04 + x_01_05
+ x_01_06 + x_01_07 + x_01_08 + x_01_09 + x_01_10 + x_01_12 + x_02_11
+ x_03_11 + x_04_11 + x_05_11 + x_06_11 + x_07_11 + x_08_11 + x_09_11
+ x_10_11 + x_11_12 >= 2
Такое ограничение нужно выделить для каждого подграфа. На первом этапе в этом примере их должно быть четыре.
Это очень похоже на то, как если бы мы взяли липкую ленту и перемотали граф в тех местах, где он пытается развалиться. Такой простой подход позволяет, нет, не в разы, а на порядки уменьшить число вычислений (для больших матриц).
Внимание! Суть метода:
Следующее ограничение мне сложно выразить математически, но звучать оно должно следующим образом:
Для каждого подмножества, на которые распадалось множество V на предыдущих этапах, должно существовать по крайней мере две связи, соединяющие его с другими подмножествами соответствующего этапа.
Процесс поиска повторять до тех пор, пока не останется только одно множество, оно и будет искомым.
По моим выкладкам получается, что даже в самом худшем случае алгоритму на симметричном полносвязном графе необходимо не больше чем шагов:
и не больше чем дополнительных ограничений:
Число 6 в сумме не случайно, на графах с числом вершин пять и менее дополнительные ограничения совсем не нужны, так как в таком случае не могут появится подциклы. (Для образования подцикла при текущей формулировке задачи требуется минимум три вершины.)
Так для n = 100, нужно не более 523, дополнительных ограничений в худшем случае, а для n=1000, 19613, что вполне не плохо.
Исходя из того, что исходные условие задаются симметричной матрицей смежности, мы можем для экономии памяти отбросить верхнюю или нижнюю её часть, а также главную диагональ.
Уменьшим, таким образом, число переменных для солвера до:
Почему же так важно уменьшить число переменных и дополнительных ограничений? Если под капотом у решателя находится симплекс метод, или его аналог, то формируется матрица:
где N – число переменных участвующих в расчёте, а M – количество ограничений.
А каждое дополнительное неравенство увеличивает исходную матрицу на 1 по обеим осям.
Чем больше матрица тем медленнее расчёты, очевидно.
(PuLP отбрасывает симметричные переменные автоматически)
Пошаговое решение граф n = 200. Осторожно много изображений!
Алгоритм так же решает здачу коммивояжёра на максимум, хоть и не столь хорошо как на минимум.
Имея решение замкнутой задачи не столь сложно немного изменить условие и сделать не замкнутый вариант.
Не убедил, что решение торт? Тестируйте!
TSP PuLP
#------------------------------------------------------------------------------------
# TSP_ILP - traveling salesman problem integer linear programming (ILP) method
# Решатель задачи комивояжора методом целочисленного линейного программирования
# Rebuilder 21.01.2023 https://habr.com/ru/publication/edit/711708/
#------------------------------------------------------------------------------------
import numpy as np
import random
import matplotlib.pyplot as plt
import networkx as nx
import pulp as pl
from datetime import datetime
#------------------------------------------------------------------------------------
def distance(x1, y1, x2, y2):
return ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5
#------------------------------------------------------------------------------------
print('Подключенные решатели:', pl.listSolvers(onlyAvailable=True))
# Количество вершин графа
n = 80
# random.seed(2)
points = [(round(random.randint(0, 1500)), round(random.randint(0, 1000))) for i in range(n)]
init_graph = nx.complete_graph(n)
for i, j in init_graph.edges():
if j > i:
init_graph[i][j]['weight'] = round(distance(points[i][0], points[i][1], points[j][0], points[j][1]))
# Развернуть граф как матрицу смежности
# adjacency_matrix = nx.adjacency_matrix(init_graph).todense()
# print(adjacency_matrix)
start_time = datetime.now()
# Создаём модель для решения
model = pl.LpProblem(name="tsp", sense=pl.LpMinimize)
#model = pl.LpProblem(name="tsp", sense=pl.LpMaximize)
# Подключаем решатель
solver = pl.PULP_CBC_CMD(msg=False) # Открытый, по умолчанию
# solver = pl.SCIP_CMD(msg=False) # Открытый, медленнее CBC
# solver = pl.CPLEX_PY(msg=False) # 45 Быстрый коммерческий (1000 Переменных)
# solver = pl.XPRESS_PY(msg=False) # ~98 Быстрый коммерческий (5000 Переменных + ограничений)
# solver = pl.GUROBI(msg=False) # 63 Быстрый коммерческий (2000 Переменных)
# Объявляем переменные
x = [pl.LpVariable(name=f'x_{i:03}_{j:03}', cat='Binary') for i in range(n) for j in range(n)]
# Задаём начальные ограничения
for i in range(n):
model += pl.lpSum([x[i * n + j if i < j else j * n + i] for j in range(n) if i != j]) == 2
# Устанавливаем целевую функцию
model += pl.lpSum([init_graph[i][j]['weight'] * x[i * n + j if i < j else j * n + i] for i in range(n) for j in range(n) if i < j])
step = 0
while True:
# Ищем решение
status = model.solve(solver)
step += 1
graph_resilt = nx.Graph()
a = 0
for i, v in enumerate(model.variables()):
# Приводим float от солвера к int
int_val = round(v.value())
if int_val > 0:
temp_name = v.name.split('_')
ii, jj = int(temp_name[1]), int(temp_name[2])
graph_resilt.add_edge(ii, jj)
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
print('Step =', step, 'Sets =', len(result_sets))
if qty_sets == 1:
break
if qty_sets == 2:
model += pl.lpSum([x[ii * n + jj if ii<jj else jj * n + ii] for i, vi in enumerate(result_sets) for j, vj in enumerate(result_sets) if i < j for ii in vi for jj in vj]) >= 2
else:
for i, vi in enumerate(result_sets):
model += pl.lpSum([x[ii * n + jj if ii<jj else jj * n + ii] for j, vj in enumerate(result_sets) if i != j for ii in vi for jj in vj]) >= 2
date_diff = (datetime.now() - start_time).total_seconds()
print('Time =', date_diff)
min_cycle = nx.find_cycle(graph_resilt, 0)
min_path = [i[0] for i in min_cycle]
# print(model)
print('Length =', round(model.objective.value()), 'steps =', step, 'path =', min_path)
plt.figure(figsize=(8, 6))
plt.axis ("equal")
nx.draw(init_graph, points, width=0, edge_color="#C0C0C0", with_labels=True, node_size=0, font_size=10)
nx.draw(graph_resilt, points, width=2, edge_color="red", style="-", with_labels=False, node_size=0, alpha=1)
plt.show()
TSP CVXPY
#------------------------------------------------------------------------------------
# TSP_ILP - traveling salesman problem integer linear programming (ILP) method
# Решатель задачи комивояжора методом целочисленного линейного программирования
# Rebuilder 21.01.2023 https://habr.com/ru/publication/edit/711708/
#------------------------------------------------------------------------------------
import numpy as np
import matplotlib.pyplot as plt
import cvxpy as cp
import random
import networkx as nx
from datetime import datetime
#------------------------------------------------------------------------------------
def distance(x1, y1, x2, y2):
# Расстояние между двумя точками на плоскости
return ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5
#------------------------------------------------------------------------------------
def to_flat_index(a, b):
# Возвращет плоский индекс для квадратной симметричной матрицы
summ = 0
if a >= b:
for i in range(1, a):
summ += i
return summ + b
else:
for i in range(1, b):
summ += i
return summ + a
#------------------------------------------------------------------------------------
def to_matrix_indexes(a):
# Возвращает кортеж индексов для квадратной симметричной матрицы по плоскому индексу
i = 0
while (a > i):
i += 1
a -= i
return (i+1, a)
#------------------------------------------------------------------------------------
def matrix_count(a):
# Возвращет размер плоского списка для квадратной симметричной матрицы
return((a * a - a) // 2)
#------------------------------------------------------------------------------------
print('Подключенные решатели:', cp.installed_solvers())
# Количество вершин графа
n = 80
# random.seed(3)
points = [(round(random.randint(0, 1500)), round(random.randint(0, 1000))) for i in range(n)]
init_graph = nx.complete_graph(n)
for i, j in init_graph.edges():
if j > i:
init_graph[i][j]['weight'] = round(distance(points[i][0], points[i][1], points[j][0], points[j][1]))
# Развернуть граф как матрицу смежности
# adjacency_matrix = nx.adjacency_matrix(init_graph).todense()
# print(adjacency_matrix)
start_time = datetime.now()
# Количество переменных
flat_size = matrix_count(n)
c = [0] * flat_size
# Устанавливаем переменные
for i in range(flat_size):
temp = to_matrix_indexes(i)
c[i] = init_graph[temp[1]][temp[0]]['weight']
x = cp.Variable(flat_size, boolean=True)
# Ограничения
constraints = []
for i in range(n):
constraints += [cp.sum([x[to_flat_index(i, j)] for j in range(n) if i != j]) == 2]
# Целевая функция
objective = cp.Minimize(cp.sum(c @ x))
step = 0
while True:
# Определяем задачу
prob = cp.Problem(objective, constraints)
# Ищем решение
prob.solve(solver = cp.GLPK_MI)
# prob.solve(solver=cp.SCIPY)
# prob.solve(solver = cp.CBC)
# prob.solve(solver = cp.XPRESS)
# print("status:", prob.status)
step += 1
graph_resilt = nx.Graph()
for i in np.argwhere(x.value == True):
temp = to_matrix_indexes(round(i[0]))
graph_resilt.add_edge(temp[1], temp[0])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
print('Step =', step, 'Sets =', len(result_sets))
# Решение найдено
if qty_sets == 1:
break
# Добавочные ограничения
if qty_sets == 2:
constraints += [cp.sum([x[to_flat_index(i, j)] for i in result_sets[0] for j in result_sets[1]]) >= 2]
else:
for i, vi in enumerate(result_sets):
constraints += [cp.sum([x[to_flat_index(ii, jj)] for j, vj in enumerate(result_sets) if i != j for ii in vi for jj in vj]) >= 2]
date_diff = (datetime.now() - start_time).total_seconds()
print('Time =', date_diff)
min_cycle = nx.find_cycle(graph_resilt, 0)
min_path = [i[0] for i in min_cycle]
print('Length =', round(prob.value), 'steps =', step, 'path =', min_path)
plt.figure(figsize=(8, 6))
plt.axis ("equal")
nx.draw(init_graph, points, width=0, edge_color="#C0C0C0", with_labels=True, node_size=0, font_size=10)
nx.draw(graph_resilt, points, width=2, edge_color="red", style="-", with_labels=False, node_size=0, alpha=1)
plt.show()
Так же предлагаю оценить собственную реализацию с использованием открытой библиотеки SciPy. В SciPy 1.10.0.. (scipy.optimize.linprog) Появился замечательный целочисленный решатель.
scipy.optimize.linprog(c, A_ub=None, b_ub=None, A_eq=None, b_eq=None, bounds=None, method='highs', callback=None, options=None, x0=None, integrality=None)
Задать его можно, скомбинировав опции: [method='highs', integrality=1].
Код выглядит сложноее. Такая задача со звёздочкой. Но и работает быстрее, чем варианты выше.
TSP SciPy
#------------------------------------------------------------------------------------
# TSP_ILP - traveling salesman problem integer linear programming (ILP) method
# Решатель задачи комивояжора методом целочисленного линейного программирования
# Rebuilder 21.01.2023 https://habr.com/ru/publication/edit/711708/
#------------------------------------------------------------------------------------
import numpy as np
from scipy import optimize
import random
import matplotlib.pyplot as plt
import networkx as nx
from datetime import datetime
#------------------------------------------------------------------------------------
def distance(x1, y1, x2, y2):
# Расстояние между двумя точками на плоскости
return ((x2 - x1) ** 2 + (y2 - y1) ** 2) ** 0.5
#------------------------------------------------------------------------------------
def to_flat_index(a, b):
# Возвращет плоский индекс для квадратной симметричной матрицы
summ = 0
if a >= b:
for i in range(1, a):
summ += i
return summ + b
else:
for i in range(1, b):
summ += i
return summ + a
#------------------------------------------------------------------------------------
def to_matrix_indexes(a):
# Возвращает кортеж индексов для квадратной симметричной матрицы по плоскому индексу
i = 0
while (a > i):
i += 1
a -= i
return (i+1, a)
#------------------------------------------------------------------------------------
def matrix_count(a):
# Возвращет размер плоского списка для квадратной симметричной матрицы
return((a * a - a) // 2)
#------------------------------------------------------------------------------------
def TSP_ILP_cycle(graph, direction=1):
# Находит минимальный (максимальный) путь в графе методом целочисленного линейного программирования
# direction - Направление оптимизации: 1 Минимизировать путь, -1 Максимизировать путь
print('Run', end='')
# Количество вершин графа
n = len(graph.nodes())
# Количество переменных
flat_n = matrix_count(n)
# Целевая функция
objective_function = np.array([direction * graph[i][j]['weight'] for i in range(n) for j in range(n) if i > j], dtype=np.int32)
# Ограничение на диаппазон переменных
bounds = (0, 1)
# Объявление массивов ограничений
equality_a = np.zeros((n, flat_n), dtype=np.int8)
equality_b = np.array([2] * n, dtype=np.int8)
inequalities_a = np.zeros((0, flat_n), dtype=np.int8)
inequalities_b = np.zeros(0, dtype=np.int8)
# Начальные ограничения модели
for i in range(n):
for j in range(n):
if i != j:
equality_a[i][to_flat_index(i, j)] = 1
# Ищем опорное решение
res = optimize.linprog(c=objective_function, bounds=bounds, A_eq=equality_a, b_eq=equality_b, method='highs', integrality=1)
step = 1
constraints = 0
while True:
print('.', end='')
# Сохраняем результат в виде графа
graph_resilt = nx.Graph()
for i, val in enumerate(res.x):
# Приводим float от солвера к int
int_val = round(val)
if int_val == 1:
temp = to_matrix_indexes(i)
graph_resilt.add_edge(temp[1], temp[0])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
# # Рисуем граф
# fig = plt.figure(figsize=(6.8, 4), edgecolor='black', linewidth=0)
# plt.title(f'Size: {n}; Length: {direction * round(res.fun)}; Steps: {step}; Constraints: {constraints}', fontsize=10)
# plt.axis("equal")
# nx.draw(graph, points, width=0, with_labels=True, node_size=0, font_size=6, font_color="black", node_shape='o')
# nx.draw(graph_resilt, points, width=1, edge_color="red", style="-", with_labels=False, node_size=0)
# Решение найдено если в графе только одно множество
if qty_sets == 1:
break
# Вводим дополнительные ограничения модели
if qty_sets == 2:
temp = np.zeros((1, flat_n), dtype=np.int8)
for i in result_sets[0]:
for j in result_sets[1]:
temp[0][to_flat_index(i, j)] = -1
inequalities_a = np.append(inequalities_a, temp, axis=0)
inequalities_b = np.append(inequalities_b, np.array([-2], dtype=np.int8))
constraints += 1
elif qty_sets > 2:
temp = np.zeros((qty_sets, flat_n), dtype=np.int8)
for i, vi in enumerate(result_sets):
for j, vj in enumerate(result_sets):
if i != j:
for ii in vi:
for jj in vj:
temp[i][to_flat_index(ii, jj)] = -1
inequalities_a = np.append(inequalities_a, temp, axis=0)
inequalities_b = np.append(inequalities_b, np.array([-2] * qty_sets, dtype=np.int8))
constraints += qty_sets
# Ищем уточненное решение
res = optimize.linprog(c=objective_function, bounds=bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1)
step += 1
print()
return {'length' : direction * round(res.fun), 'constraints' : constraints, 'steps' : step, 'path' : [i[0] for i in nx.find_cycle(graph_resilt, 0)], 'graph' : graph_resilt}
#-----------------------------------------------------------------------------------
def TSP_ILP_path(graph, end_points=(), direction=1):
# Находит минимальный (максимальный) путь в графе методом целочисленного линейного программирования
# direction - Направление оптимизации: 1 Минимизировать путь, -1 Максимизировать путь
def get_end_points(li, n):
# Возвращает конечные точки пути
path = [0] * n
for i in li:
for j in i:
path[j] += 1
return [i for i, val in enumerate(path) if val == 1]
print('Run', end='')
# Количество вершин графа
n = len(graph.nodes())
# Количество переменных
flat_n = matrix_count(n)
# Определяем число заданных граничных вершин
end_points_n = 0
if len(end_points) > 0 and end_points[0] >= 0 and end_points[0] < n:
end_points_n += 1
if len(end_points) > 1 and end_points[1] >= 0 and end_points[1] < n and end_points[0] != end_points[1]:
end_points_n += 1
# Целевая функция
objective_function = np.array([direction * graph[i][j]['weight'] for i in range(n) for j in range(n) if i > j], dtype=np.int32)
# Ограничение на диаппазон переменных
bounds = (0, 1)
# Объявление массивов ограничений
equality_a = np.ones((1, flat_n), dtype=np.int8)
equality_b = np.zeros(1, dtype=np.int32)
inequalities_a = np.zeros(((n - end_points_n) * 2, flat_n), dtype=np.int8)
inequalities_b = np.zeros((n - end_points_n) * 2, dtype=np.int8)
# Начальные ограничения модели
equality_b[0] = n - 1
k = 0
for i in range(n):
# Если были заданы конечные точки, то исключаем их из неравенств для уменьшения условий
if (end_points_n > 0) and (end_points[0] == i) or (end_points_n > 1) and (end_points[1] == i):
continue
for j in range(n):
if i != j:
flat_index = to_flat_index(i, j)
inequalities_a[k][flat_index] = -1
inequalities_a[k + 1][flat_index] = 1
inequalities_b[k] = -1
inequalities_b[k + 1] = 2
k += 2
# Если были заданы конечные точки, то добавляем условиие на равенство единице в вершинах таких точек
if end_points_n > 0:
equality_a = np.append(equality_a, np.zeros((1, flat_n), dtype=np.int8), axis=0)
equality_b = np.append(equality_b, np.ones(1, dtype=np.int8))
for i in range(n):
if i != end_points[0]:
equality_a[1][to_flat_index(i, end_points[0])] = 1
if end_points_n > 1:
equality_a = np.append(equality_a, np.zeros((1, flat_n), dtype=np.int8), axis=0)
equality_b = np.append(equality_b, np.ones(1, dtype=np.int8))
for i in range(n):
if i != end_points[1]:
equality_a[2][to_flat_index(i, end_points[1])] = 1
# Ищем опорное решение
res = optimize.linprog(c=objective_function, bounds=bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1)
step = 1
constraints = 0
while True:
print('.', end='')
# Сохраняем результат в виде графа
graph_resilt = nx.Graph()
for i, val in enumerate(res.x):
# Приводим float от солвера к int
int_val = round(val)
if int_val == 1:
temp = to_matrix_indexes(i)
graph_resilt.add_edge(temp[1], temp[0])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
# Решение найдено если в графе только одно множество
if qty_sets == 1:
break
# Вводим дополнительные ограничения модели
if qty_sets == 2:
temp = np.zeros((1, flat_n), dtype=np.int8)
for i in result_sets[0]:
for j in result_sets[1]:
temp[0][to_flat_index(i, j)] = -1
inequalities_a = np.append(inequalities_a, temp, axis=0)
inequalities_b = np.append(inequalities_b, np.array([-2], dtype=np.int8))
constraints += 1
elif qty_sets > 2:
temp = np.zeros((qty_sets, flat_n), dtype=np.int8)
for i, vi in enumerate(result_sets):
for j, vj in enumerate(result_sets):
if i != j:
for ii in vi:
for jj in vj:
temp[i][to_flat_index(ii, jj)] = -1
inequalities_a = np.append(inequalities_a, temp, axis=0)
inequalities_b = np.append(inequalities_b, np.array([-2] * qty_sets, dtype=np.int8))
constraints += qty_sets
# Ищем уточненное решение
res = optimize.linprog(c=objective_function, bounds=bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1)
step += 1
# Отображаем минимальный путь
ep = get_end_points(graph_resilt.edges, n)
path = list(nx.all_simple_paths(graph_resilt, ep[0], ep[1]))[0]
print()
return {'length' : direction * round(res.fun), 'constraints' : constraints, 'steps' : step, 'path' : path, 'graph' : graph_resilt}
#-----------------------------------------------------------------------------------
# Количество вершин графа
n = 80
# random.seed(4)
points = [(round(random.randint(0, 15000)), round(random.randint(0, 10000))) for i in range(n)]
init_graph = nx.complete_graph(n)
for i, j in init_graph.edges():
if j > i:
init_graph[i][j]['weight'] = round(distance(points[i][0], points[i][1], points[j][0], points[j][1]))
# Развернуть граф как матрицу смежности
# adjacency_matrix = nx.adjacency_matrix(init_graph).todense()
# print(adjacency_matrix)
start_time = datetime.now()
# Запускаем расчёт замкнутой TSP
res = TSP_ILP_cycle(init_graph, direction=1)
# Запускаем расчёт незамкнутой TSP
# res = TSP_ILP_path(init_graph, end_points=[], direction=1)
date_diff = (datetime.now() - start_time).total_seconds()
print('Time =', date_diff)
print('Path =', res['path'])
# Рисуем граф
fig = plt.figure(figsize=(6.8, 4), edgecolor='black', linewidth=0)
plt.title(f'Size: {n}; Length: {res["length"]}; Steps: {res["steps"]}; Constraints: {res["constraints"]}', fontsize=10)
plt.axis("equal")
nx.draw(init_graph, points, width=0, with_labels=True, node_size=0, font_size=6, font_color="black", node_shape='o')
nx.draw(res["graph"], points, width=1, edge_color="red", style="-", with_labels=False, node_size=0)
Теперь о грустном. Вся сложность вычислений ложится на решатель линейных уравнений. От его качества зависит практически всё. Как я понимаю, в среднем решение находится как O(n^3) от числа вершин. К сожалению, есть и плохие случаи, когда решение падает до O(2^n). Целочисленное программирование является NP-трудной задачей. Следовательно задача коммивояжёра остаётся NP-трудной, несмотря на то, что мы немного снизили асимптоматику задачи.
Эвристические методы (генетический и муравьиный алгоритмы) конечно обходят текущую реализацию алгоритма по скорости нахождения решения, но зато грешат в точности.
Что даёт нам точное решение проблемы коммивояжёра? Как мне видится, кроме очевидного решения транспортной задачи, наибольший выигрыш может получить схемотехника / материалообработка.
На печатной плате нужно просверлить монтажные отверстия (произвести пайку / прочее действие), в какой последовательности нужно обходить точки контакта, чтобы максимально сократить время на перемещение рабочего инструмента? Задача коммивояжёра в чистом виде.
Полагаю, что сам подход полученный в ходе решения данной задачи может помочь найти ключ к решению других задач класса NP.
P.S. Создать эту работу было NP-трудно, надеюсь она вам понравилась.
Комментарии (31)
KongEnGe
21.01.2023 11:45+3Очень не хватает формального описания подхода для понимания сути происходящего.
rebuilder Автор
21.01.2023 12:24+1Очень не хватает формального описания подхода для понимания сути происходящего.
Теоретическое обоснование разжевано десятилетия назад. Читать не математику это будет не интересно, а тот к то в теме поймёт опираясь на свой опыт работы с TSP.
Всё что я сделал, это применил новый подход к выбору ограничивающих условий.
KongEnGe
21.01.2023 12:46+2Тем более, если теоретически обосновано, почему не дать ссылку на готовое? А то я к TSP небезразличен, но вот реконструировать подход по косвенным признакам уже лень :)
То, что я вижу, никак не говорит о том, что в варианте задачи дискретной оптимизации есть уход от NP-полноты. Главное матрицу расстояний попротивнее подогнать на обсчет.
rebuilder Автор
21.01.2023 13:02Понимаете коллега, где вы тут увидели дискретную оптимизацию? В ней я разочаровался еще, когда делал ветви и границы. Как раз линейное целочисленное программирование, основанное на симплекс методе и алгоритме Гомори, как раз и даёт эти чудесные результаты.
KongEnGe
21.01.2023 13:08+1У Вас после жирной "Сути метода" пошла формула как раз из раздела "Дискретная оптимизация":
https://ru.wikipedia.org/wiki/Задача_коммивояжёра#Формулировка_в_виде_задачи_дискретной_оптимизацииПоэтому и хотелось увидеть, где случился ранее незамеченный поворот в прорывную методику.
rebuilder Автор
21.01.2023 13:26О! Я так рад, что вы заинтересовались и начали копать. Вы абсолютно правы, за основу была выбрана именно эта формулировка. Но если посмотреть работы Данцига, Фалкерсона и Джонсона, а так же Миллера, Такера и Землина (ссылки на эти работы есть там же в википедии), то можно прийти к следующему выводу:
Эти авторы пытались засунуть все ограничения в систему уравнений сразу, а так их получается очень много и решатель захлёбывается. Я же предлагаю добавлять только те ограничения, которые нужны по факту в конкретной задаче.
Надеюсь так яснее. Я думал код приведённый выше сам за себя скажет.
rebuilder Автор
21.01.2023 15:12Тут я не совсем корректно выразился. Выше названные исследователи, конечно же, перебирали комбинации. Если выразить эти комбинации в виде линейных неравенств и отдать их целочисленному решателю, то будет слишком много уравнений.
KongEnGe
22.01.2023 11:45Мне просто все происходящее напомнило одну давнюю дискуссию про "коммивояжера за полином" (https://forum.ixbt.com/topic.cgi?id=40:3025)
Там была и ссылка на наборы "нехороших" данных для задачи, которая, увы, протухла, но, вроде бы, я нашел именное ее файл в своей домашней помойке. Не хотите ли скормить в свой алгоритм?
https://drive.google.com/file/d/11PfpGmwgUCbyg9xeyEeQYpa-niWp53sL
rebuilder Автор
22.01.2023 12:22Не хотите ли скормить в свой алгоритм?
Только если вы расскажете, как это прочитать. В архиве что угодно, но только не симметричная матрица смежности.
По поводу дискуссии на форуме, как мне показалось, там речь шла про точность самой вводимой матрицы.
В тексте статьи был обозначен худший случай, требующий (n/2+1) вызовов решателя. В своих изысканиях специально искал примеры, которые являются адом для решателя, но и там он справлялся не плохо. Хоть и нелегко, но вдвое бил динамическое программирование по размеру матрицы. Если у динамического программирования сложность O(n^2*2^n), то тут я предполагаю она упирается в худшем случае в O(n*2^n) но это прикидочная величина.
KongEnGe
22.01.2023 12:33Не скажу за все, но некоторые наборы симметричные при беглом взгляде (br17). Матрицу вычитывать последовательно, значение за значением, в количестве, определяемом в DIMENSION
guyfawkes
23.01.2023 23:24Простите, но ваш код невозможно читать. Вложенные циклы, переменные с ничего не говорящими названиями. Что происходит, без отладчика совершенно неясно, особенно если человек не понимает вашего хода мысли. Общие моменты текста статьи ещё кое-как маппятся на вашу реализацию, но детали остаются неочевидны.
rebuilder Автор
24.01.2023 09:17Признаю код не самый простой, есть грешок. Самые сложные места постарался прокомментировать по ходу. Вложенные циклы можно спокойно развернуть и добавить точки отладчика чтобы посмотреть сложные места. Непонятные переменные только для циклов, для краткости, не люблю, знаете, длинные строки.
Если есть конкретные вопросы спрашивайте отвечу. Так долго варился в этом коде что не могу критически оценить его сложность.
thevlad
21.01.2023 23:36+2Дискретная оптимизация подразумевает различные методы. В данном случаи используется формулировка, как задачи целочисленного программирования. Прелесть в том, что матрицу подогнать конечно теоретически можно, но на практике это довольно сложно. Уже давал в предыдущей статье автора ссылку на SOTA решатель, основанный на тех же принципах https://en.wikipedia.org/wiki/Concorde_TSP_Solver, он может решать инстансы до десятков тысяч вершин. В то время как другие наивные подходы уже взрываются на паре десятков.
rebuilder Автор
21.01.2023 23:45Благодарю вас за наводку на солвер Concorde, он лучший. Однако после ознакомления с их описание задачи у меня сложилось чёткое убеждение, что они не используют именно целочисленное линейное программирование, скорее там обычное линейное программирование. А так же я не увидел именно точного решения. Возможно, вы мне поможете с этим разобраться?
thevlad
22.01.2023 00:23+1У них там точное решение https://www.reddit.com/r/compsci/comments/8auwm9/how_does_concorde_claim_to_be_a_tsp_solver/ В частности "Concorde is an exact algorithm based on Dantzig's cutting plane method which has obtained an optimal tour for an instance with 85900 cities." Впечатление возможно сложилось потому, что они неявно реализуют решение целочисленной формулировки через https://en.wikipedia.org/wiki/Branch_and_cut используя релаксацию. Весь дьявол в деталях реализации branc-and-cut, но я сильно глубоко не разбирался.
PS: в частности есть книжка от автора Concorde "The Traveling Salesman Problem: A Computational Study" там есть много подробностей, она лежит на libgen
tessob
23.01.2023 17:18+1Было бы интересно взглянуть на то как формулируются линейные ограничения для задачи коммивояжера. Просто из статьи это никак не следует, а сам по себе Симплекс - ни разу не магия.
abutorin
21.01.2023 12:09Из открытых библиотек есть еще GLPK, OR-Tools от Google.
rebuilder Автор
21.01.2023 12:17+1GLPK у меня не захотел решать целочисленную задачу, пришлось забраковать.
А какой решатель встроен OR-Tools, по умолчанию?
abutorin
21.01.2023 12:24А какой решатель встроен OR-Tools, по умолчанию?
По умолчанию у них свой решатель. Но можно и другие подключить. Правда другие я не пробовал. На линейних задачах встроенный быстрее чем GLPK. В свое время проводил сравнение SciPy. GLPK, OR-Tools по скорости и максимальной размерности задачи. SciPy отказался решать систему с 5000 переменных. GLPK и OR-Tools с ними справляются, последний несколько быстрее.
rebuilder Автор
21.01.2023 12:36КДПВ для статьи сформирована SciPy и для неё потребовалось 280875 переменных, нормально справился. Обратите внимание, что в последних версиях scipy.optimize.linprog завезли классные быстрые решатели. Особенно меня порадовал метод глубинной точки, для нецелочисленных задач.
abutorin
21.01.2023 13:04для неё потребовалось 280875 переменных
Тут наверное пора уточнить что такое "переменная". Я подрузумелвал задачу оптимизации матрицы в которой больше 5000 строк и столбцов, т.е. там 25 млн. ячеек.
Но спасибо за наводку, посмотрю свежую версию SciPy.
rebuilder Автор
22.01.2023 15:49Стандартный решатель OR-Tools, меня совсем не порадовал, по скорости он на последнем месте.
DonAgosto
21.01.2023 15:31На печатной плате нужно просверлить монтажные отверстия (произвести пайку / прочее действие)
Заметная польза будет, если нахождение в «вершине» дешевле, чем перемещение по «ребру». К сожалению чаще бывает наоборот.rebuilder Автор
21.01.2023 16:55Всё так, но ведь и операции могут быть разными. Например, монтажный пистолет работает достаточно шустро, прозвонка дорожек мультиметром то же весьма скоростная операция. К сожалению, не владею реальными кейсами.
tzlom
22.01.2023 12:07Примерно так работает Active Set Solver для QP проблем.
Правда в нем еще есть шаг релаксации ограничений, т.к. можно в процессе построить ограничения которые будут загонять в локальный минимум.
Кстати для AS солвера мне помогло не решать проблему полностью заново, а дорешивать её для новых условий , сложность падает с n^3 до ~n^2
Интересно, что будет если уточняющие ограничения вводить с условием >=1 ?
rebuilder Автор
22.01.2023 12:37Интересно, что будет если уточняющие ограничения вводить с условием >=1 ?
Решение всё равно найдётся но, большим числом неравенств
tzlom
22.01.2023 12:41Это теория или попробовали? Ограничение на замкнутость должно дополнять это ограничение в любом случае, разве нет?
rebuilder Автор
22.01.2023 12:47+1Попробовал, множества всё равно связываются хоть и одной связью. За счёт повторов связывания, рано или поздно стягивает всю конструкцию в монолит. Но получается медленнее чем с двумя связями.
Krasnoarmeec
Статья супер!
Хабр опять торт!
Особенно впечатляет визуализация решений при разном количестве ограничений.
Сильно жаль, что могу поставить за статью только два плюсика. ????