Задача коммивояжёра — ещё немного больше, ещё немного быстрее
Всегда выбирайте самый трудный путь — на нем вы не встретите конкурентов.
Шарль де Голль
И снова здравствуйте, уважаемые читатели Хабра. Мы продолжаем наше путешествие в мир алгоритмов поиска оптимального пути.
В прошлой работе мы уже узнали, как можно найти оптимальный путь в графе в несколько сотен вершин. В данной работе хочу более подробно остановится на сути метода, а также разобрать возможность по его ускорению на графах от тысячи элементов.
Задача коммивояжёра обычно задаётся симметричной матрицей смежности размера n x n, где n — количество вершин для графа, для которого мы хотим найти оптимальный путь.
[[ 0 106 71 81 75 140 108 146 262 174]
[106 0 162 184 92 130 36 157 270 204]
[ 71 162 0 38 88 134 151 123 223 130]
[ 81 184 38 0 121 171 179 161 257 164]
[ 75 92 88 121 0 65 69 77 196 118]
[140 130 134 171 65 0 96 33 140 84]
[108 36 151 179 69 96 0 124 235 173]
[146 157 123 161 77 33 124 0 119 52]
[262 270 223 257 196 140 235 119 0 93]
[174 204 130 164 118 84 173 52 93 0]]Для метрической задачи на плоскости такую матрицу можно легко построить, рассчитав расстояние между каждой парой точек.
def distance(point1, point2):
return math.hypot(point2[0] - point1[0], point2[1] - point1[1])Алгоритм рассмотренный далее решает не только метрическую задачу, и не обязательно в двумерном пространстве, но для удобства тестирования мы будем применять случайный граф на плоскости. Основное условие, чтобы задающая матрица была симметричной относительно главной диагонали, а сама главная диагональ игнорируется.
При таком подходе мне больше нравится задавать входное условие через полно связный неориентированный граф:
import networkx as nx
init_graph = nx.complete_graph(n)Подсчитать количество возможных рёбер для такого графа можно с помощью функции:
def matrix_count(a):
# Возвращет количество рёбер квадратной симметричной матрицы
return (a * a - a) // 2
При n = 1000 у нас получится 499500 рёбер, среди которых мы должны искать оптимальный путь.
Это очень много для задачи линейного программирования (LP), ведь каждое такое ребро это дополнительные переменная участвующая в расчётах.
Для расчёта такой задачи неявно строится симплекс таблица размером n' x m', где m' > n', n' — число участвующих переменных, m' — число неравенств. (Здесь подразумевается решение с неравенствами и симплекс метод с искусственным базисом. Хотя возможно использовать метод с естественным базисом, так можно сэкономить немного памяти на вспомогательных переменных, но принципиальной разницы нет.)
Матрица для задачи LP с таким количеством переменных занимает всю оперативную память компьютера, а это десятки Гбайт.
Мне удалось найти только два варианта разрешения ситуации с объёмом используемой памяти:
Нахождение возможности более компактного представления матрицы в памяти.
Уменьшение числа используемых переменных для уменьшения размера матрицы.
Нахождение возможности более компактного представления матрицы в памяти.
Если посмотреть на матрицу ограничений ниже, то можно увидеть, что большинство элементов равны нулю. Для хранения таких матриц можно использовать разряженные матрицы (Разрежённая матрица — матрица с преимущественно нулевыми элементами).
[[1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 1 1 1 0 0 0]
[1 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 1 1 0]
[0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 1]
[0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 1 1 1 1 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 1 0 1 1]
[0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0]
[0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0]
[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]
[0 0 1 1 0 0 1 0 0 0 1 1 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 1 1 1 0 1 0 1 0 0 1 0 0 1 1 0]]В разряженной матрице сохраняются только не нулевые элементы, как правило в виде списка списков.
К счастью в python в пакете scipy существует не только удобная библиотека для работы с разряженными матрицами, но и решатели задач линейного программирования с поддержкой таких матриц.
Существует два решателя подходящих для наших целей, они имеют небольшие нюансы по заполнению исходных данных, но в целом я не заметил весомых отличий в пользу одного из них, как по скорости работы, так и по удобству применения.
a) scipy.optimize.linprog
b) scipy.optimize.milp
Оба решателя при правильном подборе параметров используют быстрый двойной решатель целочисленных линейных уравнений с применением разряженных матриц.
res = optimize.linprog(c=objective_function, bounds=variable_bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1, options={'disp':True})
res = optimize.milp(c=objective_function, bounds=variable_bounds, constraints=optimize.LinearConstraint(matrix_a, lower_bound, upper_bound), integrality=1, options={'disp':True})Отличие только в том, что linprog разбивает условия на два блока: равенства и неравенства. Milp же использует один блок неравенств с верхним и нижним ограничением (для равенств верхнее и нижнее ограничение делается одинаковым).
Уменьшение числа используемых переменных для уменьшения размера матрицы.
Предлагаю использовать эвристику, для удаления рёбер, которые не могут участвовать в построении оптимального пути. Так можно уменьшить объём используемой памяти, а также значительно ускорить расчёт.
Давайте посмотрим на следующий произвольный граф.
Полносвязный граф, n = 100
Возьмём случайное ребро и спросим себя: может ли оно входить в оптимальный путь?
Моя биологическая нейросеть подсказывает, что данное конкретное ребро не может, но есть ли формальное описание этого правила. У меня нет точного ответа на данный вопрос, однако в ходе сотен экспериментов с поиском пути на графах можно заметить определённые закономерности. Предлагаю ниже несколько метафор объясняющих свои рассуждения.
Представим, что вершины графа — это центры гравитации — «звёзды», чем они ближе друг к другу тем более вероятно, что их общее ребро будет входить в оптимальный путь.
Каждая такая звезда может испустить строго два «луча». По одному из них мы приходим в вершину, по-другому выходим. Так как граф неориентированный нам не важно направление луча.
Если представить, что звёзды имеют не нулевой радиус, и могут отбрасывать тень на другие вершины, находящиеся за собой, то можно использовать приём очень похожий на трассировку лучей (Ray Tracing) для удаления неоптимальных рёбер.
Мы будем полагать, что по условию задачи нам не будут известны координаты вершин, а также метрика пространства, в котором развёрнут граф (известна только длина рёбер, заданная матрицей смежности). Поэтому попробуем создать очень упрощённую модель пересечения лучей.
Первой попыткой стала реализация зон контроля для каждой вершины. Радиусы подобраны так, чтобы максимизировать площадь окружностей, при этом не допустив перекрытия областей контроля.
Зоны контроля для вершин
Окружности создают вымышленную «тень». Объекты, находящиеся в геометрической тени от окружности, как правило не могут организовать связь с целевой вершиной и могут быть вычеркнуты в матрице смежности. Однако это не совсем так, в некоторых случаях возможна ситуация, когда вершина находится в тени, но ребро входит в оптимальный путь. Это связанно с тем, что для вершины всегда должны быть ровно две связи с соседними вершинами. Если удачных вариантов нет, оптимальным вариантом будет ребро, связанное с видимой вершиной видимой вершины (метафорически, луч отражённый от вершины стремится к своим соседям).
Формализуем. В полносвязном графе заданной матрицей смежности оставляем только рёбра, связанные с видимой вершиной, а также рёбра из текущей вершины в вершины, которые видны из видимой вершины. Остальные рёбра отбрасываем и не используем в расчётах, так как они не могут войти в оптимальный путь. (Понимаю, тут спорно, чёткого доказательства у меня нет.)
Такой простой эвристикой можно уменьшить число переменных для расчёта в 4–6 раз.
Как можно определить, является ли вершина видимой из текущей вершины, если мы не знаем их координат? Нужно для каждой пары вершин AB проверить не существует ли третья вершина C с радиусом r, таким, что для треугольника ABC выполнятся следующие условия: все рёбра не нулевые, AB — большая сторона, высота CH > r. Если для произвольного С, условия выполняются, то будем считать, что вершина С перекрывает AB.
Возможно ли ещё уменьшить число переменных? Оказывается, что да. Для этого нужно в мысленном эксперименте представить, что «центры масс» могут взаимно проникать друг в друга.
F (1,1)
В такой интерпретации радиус окружностей соответствует расстоянию до ближайшей вершины. Заметно более плотное заполнение пространства между вершинами, дающее более плотную тень, а следовательно, более сильное отсечение отдалённых рёбер. Однако, в такой концепции отсечения мы должны разрешить «лучу» пересекать окружность один раз, не считая его угасшим.
Если организовать процесс выделения окружностей и количество разрешённых пересечений в виде функции, то можно комбинировать большие размерности с ещё лучшими характеристиками отсечения.
F (Число минимальных окружностей для вершины, число пересечений до обрыва луча)
F (2,2)
Эмпирическим путём были выделены оптимальные варианты для такой функции с одной стороны, максимально уменьшающих число переменных симплекс метода с другой, не исключающих оптимальное решение: F (1,1), F (2,2), F (3,4), F (4,5), F (5,7).
Наилучшим образом показал себя вариант F (4,5). Лучшее соотношение времени расчётов и числа переменных, поэтому в основном и использовался для тестов. (Возможны иные варианты, не учтённые автором).
Применяя такую эвристику, в среднем для случайного набора точек нужно в 20–50 раз меньше переменных решателя линейных уравнений, что значительно ускоряет расчёты.
Радиусы окружностей — это длины самых коротких рёбер для вершины. Расcчитать их невероятно просто:
# Разворачиваем граф как матрицу смежности
adjacency_matrix = nx.adjacency_matrix(init_graph, dtype=np.float64).todense()
# Сортируем матрицу смежности по строкам
adjacency_matrix.sort(axis=1)
# Выбираем только часть матрицы
vertex_array = adjacency_matrix[:, 1:circles + 1]Несколько слов про алгоритм поиска минимального пути. В прошлой статье я кратко обрисовал, как это происходит. Тут мы остановиться на этом подробнее.
Решатель линейных уравнений предназначен для поиска оптимального решения системы оных уравнений.
Первое, что нам нужно сделать, это определить целевую функцию. В данном случае, целевая функция решателя, это список длин рёбер графа, среди которых мы ищем минимальный путь.
objective_function — коэффициенты длины каждого ребра
Следующим шагом мы должны определить в виде уравнений условие, что в каждую вершину должно вести ровно два пути. Для этого мы определим столько равенств двум, сколько у нас вершин в графе.
equality_a — коэффициенты левой части уравнений
equality_b — коэффициент правой части, тут все равны 2
Третьим шагом мы должны определить ограничение неотрицательности на все переменные. Было бы очень странно, иметь путь отрицательной длины. Так как у нас будет выполнятся условие целочисленности переменных и выполнение второго шага, то возможно два варианта условия.
а) variable_bounds = (0, 1) — возможные варианты 0, 1 (ребро либо выберется, либо нет)
б) variable_bounds = (0, None) — возможные варианты 0, 1, 2 (ребро либо выберется, либо нет, либо выберется дважды)
variable_bounds = (0, 1)
variable_bounds = (0, None)
Так как ограничение сверху на каждую переменную неявно создаёт дополнительные ограничение в решателе линейных уравнений, использование первого варианта немного замедляет решатель на каждом шаге. При выборе второго варианта, решения проходят быстрее, но шагов решателя получается немного больше, так же больше получается дополнительных ограничений. Дополнительные ограничения явно добавляются на дальнейших итерациях алгоритма.
Какую конструкцию лучше использовать не определился. В разных случаях выигрывает то одно, то другое.
Следующим шагом объявим объекты для неравенств. Пока они будут пусты, мы заполним их на следующих итерациях при вводе дополнительных ограничений.
inequalities_a — коэффициенты левой части неравенств
inequalities_b — коэффициенты правой части неравенств
У нас есть все для того чтобы получить первое решение.
res = optimize.linprog(c=objective_function, bounds=variable_bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1)
# Сохраняем результат в виде графа связей
graph_resilt = nx.Graph()
for i, val in enumerate(res.x):
if round(val) >= 1:
graph_resilt.add_edge(*array_index[i])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)Если в полученом графе будет содержатся больше одного подмножества, то нам нужно внести дополнительные ограничения в задачу и выполнить повторный поиск решения.
Связать подмножества в единое целое можно двумя эквивалентными способами:
Наложить ограничения между подмножествами, стягивающими их словно скотч.
Каждое подмножество должно соединятся с другими подмножествами по крайней мере двумя связями.
На каждое из подмножеств наложить внутреннее ограничение, заставляющее его искать недостающее звено в соседнем подмножестве.
Для каждого подмножества число связей в нём должно быть, как минимум на одно меньше количества вершин в данном подмножестве.
Так как в первом случае знак условия получится ≥, а алгоритм работает только с условиями ≤, то нужно в левой и правой части уравнения изменить знак на минус. Для второго способа это не требуется, так как знак условия неравенства остаётся ≤.
Мы в праве выбирать для каждого подмножества одно из условий (или даже оба, если угодно), даже в произвольном порядке. Однако я предлагаю выбирать то из условий, которое создаст меньше значений в разряженной матрице дополнительных ограничений. Так удаётся получить некоторый выигрыш в производительности алгоритма.
Собственно нам остаётся повторять решение до тех пор, пока в графе не останется только одно множество, которое и будет искомым минимальным путём.
На каждой итерации алгоритма в систему уравнений добавляются новые условия, из-за этого каждый новый шаг получается немного медленнее предыдущего, это тоже стоит учитывать.
По сравнению с алгоритмом, предложенным в прошлой статье, данный вариант значительно выигрывает как по скорости работы, так и по способности обрабатывать графы больших размеров.
Оцениваю среднюю сложность работы алгоритма на наборе случайных точек приблизительно как O (n**3), что проигрывает Concorde, который оценивается как O (n**2.2).
Несколько слов про Concorde, это потрясающий алгоритм. Идеи, применённые в нём, дают невероятные результаты. К сожалению, мне не удалось найти исходный код реализации. Выражаю огромную благодарность William J. Cook за его выступления, он просто невероятен.
Хотя общие моменты понятны, собственная реализация прототипа Конкорд, оказалась очень сложной и медленной, гораздо хуже алгоритма, описанного в данной работе. Думаю, что для быстрой реализации нужно суметь построить свой оптимальный решатель линейных уравнений, так как нужен шаг с применением алгоритма Гомори. Все готовые решатели, которые мне попадались не разрешают использовать промежуточные данные для внесения условий целочисленности после нахождения не целочисленного решения. Без этого каждый перезапуск решателя занимает уйму времени, что в итоге делает саму идею реализации бессмысленной. Кроме того, в готовых реализациях солверов часто применяется не алгоритм Гомори, а Branch and cut, а как мне показалось, он медленнее когда ветвлений очень много.
Данный труд представляет собой авторское переосмысление общедоступной литературы по теме поиска пути и не претендует на безоговорочную истинность. Если вам известно больше, поправьте автора, чтобы рассеять его невежество.
Спасибо за внимание!
Код
Исходный код распространяется под лицензией GNU GPL
#------------------------------------------------------------------------------------
# TSP - traveling salesman problem integer linear programming method
# Решатель задачи комивояжёра методом целочисленного линейного программирования
# Rebuilder 09.06.2023 https://habr.com/ru/publication/edit/740984/
#------------------------------------------------------------------------------------
import numpy as np
import random
import matplotlib.pyplot as plt
from scipy import optimize
from scipy import sparse
import networkx as nx
from datetime import datetime
import math
#------------------------------------------------------------------------------------
def distance(point1, point2):
# Возвращет дистанцию между двумя точками на плоскости
return math.hypot(point2[0] - point1[0], point2[1] - point1[1])
#------------------------------------------------------------------------------------
def matrix_count(a):
# Возвращет размер плоского списка для квадратной симметричной матрицы
return((a * a - a) // 2)
#------------------------------------------------------------------------------------
def triangle_intersection(a, b, c, vertex_array, val):
# Нужно ли убрать ребро
if (b > a) or (c > a) or (a < 0.000001):
return val
cross_size = len(vertex_array)
key = [True] * cross_size
for i in range(cross_size):
if (b == vertex_array[i]):
key[i] = False
p = 0.5 * (a + b + c)
s = p * (p - a) * (p - b) * (p - c)
r = 2 * np.sqrt(abs(s)) / a
for i in range(cross_size):
if key[i] and vertex_array[i] >= r:
val += 1
return val
#------------------------------------------------------------------------------------
def algorithm_evaluation(array_key, result_sets, k):
# Возвращет оценку какой алгоритм создаст меньше единиц
set1 = set()
set2 = set()
for i in result_sets[k]:
for j in result_sets[k]:
if i > j and (i, j) in array_key:
set1.add(array_key[(i, j)])
for l, val_l in enumerate(result_sets):
if l != k:
for i in result_sets[k]:
for j in val_l:
key = (i, j) if i > j else (j, i)
if (key in array_key):
set2.add(array_key[key])
if len(set1) <= len(set2):
return True, set1
else:
return False, set2
# ====================================================================================================================
def TSP_linprog(init_graph, circles, cross):
# Находит минимальный путь в графе методом целочисленного линейного программирования
print('Run TSP linprog', end='')
n = len(init_graph)
# Количество переменных солвера
flat_n = matrix_count(n)
# Разворачиваем граф как матрицу смежности
adjacency_matrix = nx.adjacency_matrix(init_graph, dtype=np.float64).todense()
# Сортируем матрицу смежности по строкам
adjacency_matrix.sort(axis=1)
# Выбираем только часть матрицы
vertex_array = adjacency_matrix[:,1:circles + 1]
del adjacency_matrix
graph_of_neighbors = nx.Graph(())
for i in range(n):
for j in range(n):
if i > j:
value = 0
for k in range(n):
if (k != i) and (k != j):
value = triangle_intersection(init_graph[i][j]['weight'], init_graph[i][k]['weight'], init_graph[j][k]['weight'], vertex_array[k], value)
if value > cross:
break
if value <= cross:
graph_of_neighbors.add_edge(i, j)
del vertex_array
graph_united = nx.Graph(())
for i, j in list(graph_of_neighbors.edges()):
for k in list(graph_of_neighbors.adj[i]):
if j != k and not graph_united.has_edge(j, k):
graph_united.add_edge(j, k)
for k in list(graph_of_neighbors.adj[j]):
if i != k and not graph_united.has_edge(i, k):
graph_united.add_edge(i, k)
for i, j in list(graph_of_neighbors.edges()):
if not graph_united.has_edge(i, j):
graph_united.add_edge(i, j)
del graph_of_neighbors
# Ограничение на диаппазон переменных
variable_bounds = (0, 1)
# variable_bounds = (0, None)
# Заполняем основные элементы
edges_list = list(graph_united.edges())
n_var = len(edges_list)
array_index = {}
array_key = {}
objective_function = np.empty(n_var, dtype=np.float64)
# Объявление массивов ограничений
equality_a = sparse.lil_matrix((n, n_var), dtype=np.int8)
equality_b = [2] * n
inequalities_a = sparse.lil_matrix((0, n_var), dtype=np.int8)
inequalities_b = []
for i, val in enumerate(edges_list):
key = (val[0], val[1]) if val[0] > val[1] else (val[1], val[0])
array_key[key] = i
array_index[i] = key
objective_function[i] = init_graph[key[0]][key[1]]['weight']
equality_a[key[0], i] = 1
equality_a[key[1], i] = 1
del graph_united
constraints = 0
step = 0
while True:
# Ищем решение
res = optimize.linprog(c=objective_function, bounds=variable_bounds, A_eq=equality_a, b_eq=equality_b, A_ub=inequalities_a, b_ub=inequalities_b, method='highs', integrality=1, options={'disp':True})
step += 1
print('.', end='')
# Сохраняем результат в виде графа связей
graph_resilt = nx.Graph()
for i, val in enumerate(res.x):
if round(val) >= 1:
graph_resilt.add_edge(*array_index[i])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
# Решение найдено если в графе только одно множество
if qty_sets == 1:
break
# Если только два множества соединяем их напрямую
if qty_sets == 2:
inequalities_a.resize((constraints + 1, n_var))
inequalities_b += [-2]
for i in result_sets[0]:
for j in result_sets[1]:
key = (i, j) if i > j else (j, i)
if key in array_key:
inequalities_a[constraints, array_key[key]] = -1
constraints += 1
continue
# Вводим столько ограничений, сколько у нас получилось множеств
inequalities_a.resize((constraints + qty_sets, n_var))
for i, val_i in enumerate(result_sets):
# Выбираем тот алгоритм, который добавит меньше значений в разряженную матрицу
value, add_set = algorithm_evaluation(array_key, result_sets, i)
if value:
inequalities_b += [len(val_i) - 1]
for j in add_set:
inequalities_a[constraints + i, j] = 1
else:
inequalities_b += [-2]
for j in add_set:
inequalities_a[constraints + i, j] = -1
constraints += qty_sets
print()
return {'length' : res.fun, 'graph' : graph_resilt, 'constraints' : constraints, 'steps' : step}
# ====================================================================================================================
def TSP_milp(init_graph, circles, cross):
# Находит минимальный путь в графе методом целочисленного линейного программирования
print('Run TSP milp', end='')
n = len(init_graph)
# Количество переменных солвера
flat_n = matrix_count(n)
# Разворачиваем граф как матрицу смежности
adjacency_matrix = nx.adjacency_matrix(init_graph, dtype=np.float64).todense()
# Сортируем матрицу смежности по строкам
adjacency_matrix.sort(axis=1)
# Выбираем только часть матрицы
vertex_array = adjacency_matrix[:,1:circles + 1]
del adjacency_matrix
graph_of_neighbors = nx.Graph(())
for i in range(n):
for j in range(n):
if i > j:
value = 0
for k in range(n):
if (k != i) and (k != j):
value = triangle_intersection(init_graph[i][j]['weight'], init_graph[i][k]['weight'], init_graph[j][k]['weight'], vertex_array[k], value)
if value > cross:
break
if value <= cross:
graph_of_neighbors.add_edge(i, j)
del vertex_array
graph_united = nx.Graph(())
for i, j in list(graph_of_neighbors.edges()):
for k in list(graph_of_neighbors.adj[i]):
if j != k and not graph_united.has_edge(j, k):
graph_united.add_edge(j, k)
for k in list(graph_of_neighbors.adj[j]):
if i != k and not graph_united.has_edge(i, k):
graph_united.add_edge(i, k)
for i, j in list(graph_of_neighbors.edges()):
if not graph_united.has_edge(i, j):
graph_united.add_edge(i, j)
del graph_of_neighbors
# Заполняем основные элементы
edges_list = list(graph_united.edges())
n_var = len(edges_list)
array_index = {}
array_key = {}
matrix_a = sparse.lil_matrix((n, n_var), dtype=np.int8)
objective_function = np.empty(n_var, dtype=np.float64)
for i, val in enumerate(edges_list):
key = (val[0], val[1]) if val[0] > val[1] else (val[1], val[0])
array_key[key] = i
array_index[i] = key
objective_function[i] = init_graph[key[0]][key[1]]['weight']
matrix_a[key[0], i] = 1
matrix_a[key[1], i] = 1
del graph_united
# Ограничение на диаппазон переменных
variable_bounds = (0, 1)
# variable_bounds = (0, np.inf)
# Начальные ограничения модели
lower_bound = [2] * n
upper_bound = [2] * n
constraints = 0
step = 0
while True:
# Ищем решение
res = optimize.milp(c=objective_function, constraints=optimize.LinearConstraint(matrix_a, lower_bound, upper_bound), bounds=variable_bounds, integrality=1, options={'disp':True})
step += 1
print('.', end='')
# Сохраняем результат в виде графа связей
graph_resilt = nx.Graph()
for i, val in enumerate(res.x):
if round(val) >= 1:
graph_resilt.add_edge(*array_index[i])
# Разбиваем граф результата на подмножества
result_sets = list(nx.connected_components(graph_resilt))
qty_sets = len(result_sets)
# Решение найдено если в графе только одно множество
if qty_sets == 1:
break
# Если только два множества соединяем их напрямую
if qty_sets == 2:
matrix_a.resize((n + constraints + 1, n_var))
lower_bound += [2]
upper_bound += [np.inf]
for i in result_sets[0]:
for j in result_sets[1]:
key = (i, j) if i > j else (j, i)
if key in array_key:
matrix_a[n + constraints, array_key[key]] = 1
constraints += 1
continue
# Вводим столько ограничений, сколько у нас получилось множеств
matrix_a.resize((n + constraints + qty_sets, n_var))
for i, val_i in enumerate(result_sets):
# Выбираем тот алгоритм, который добавит меньше значений в разряженную матрицу
value, add_set = algorithm_evaluation(array_key, result_sets, i)
if value:
lower_bound += [0]
upper_bound += [len(val_i) - 1]
else:
lower_bound += [2]
upper_bound += [np.inf]
for j in add_set:
matrix_a[n + constraints + i, j] = 1
constraints += qty_sets
print()
return {'length' : res.fun, 'graph' : graph_resilt, 'constraints' : constraints, 'steps' : step}
# ====================================================================================================================
n = 200
# random.seed(1)
points = [(round(random.randint(0, 300)), round(random.randint(0, 200))) 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'] = distance(points[i], points[j])
plt.figure(figsize=(12, 9))
plt.axis("equal")
start_time = datetime.now()
res = TSP_linprog(init_graph, 4, 5)
# res = TSP_milp(init_graph, 4, 5)
date_diff = (datetime.now() - start_time).total_seconds()
print('time =', date_diff)
nx.draw(init_graph, points, width=0, with_labels=True, node_size=0, font_size=6, font_color="black")
nx.draw(res["graph"], points, width=1, edge_color="red", style="-", with_labels=False, node_size=0)
plt.show()
# ====================================================================================================================Пошаговое решение графа, n = 1200 (Осторожно трафик)
