Все пути одинаковы: они ведут в никуда. Но у одних есть сердце, а у других — нет. Один путь дает тебе силы, другой — уничтожает тебя.

- Карлос Кастанеда

Здравствуйте уважаемые дамы и господа, а также не бинарные личности. Хорошей эпохи.

Мы с вами уже пробовали решать точно задачу коммивояжёра методом динамического программирования и методом ветвей и границ, результат не плох, но слабоват.

В данной статье постараюсь показать, что точное решение ближе, чем принято считать.

Мы будем использовать метод целочисленного программирования, который является частным случаем линейного программирования, который в свою очередь является подклассом математического программирования.

Математическое программирование - это область математики, разрабатывающая теорию, численные методы решения многомерных задач оптимизации с ограничениями.

Ограничения могут быть:

  • дискретными (дискретное программирование или комбинаторная оптимизация (метод динамического программирования, ветвей и границ можно отнести сюда))

  • линейными (целочисленными или непрерывными)

  • нелинейными

Но и сама целевая функция может быть:

  • линейной (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.

  1. PuLP - даже из названия становится понятно, что это библиотека, которая предназначена для решения задач линейного программирования. Среди иных альтернатив она выделяется минимализмом, проста для освоения новичками, отличается удобством использования и возможностью подключать внешние решатели.

  2. CVXPY - это встроенный в Python язык моделирования с открытым исходным кодом для задач выпуклой оптимизации. В отличие от PuLP, может решать не только линейные задачи. Содержит более богатый функционал и приличный набор подключаемых решателей. Однако, как мне показалось CVXPY справляется с однотипными задачами медленнее, чем PuLP на тех же внешних решателях (по крайней мере, MIP) и более сложен для изучения. Если вам не требуется для работы выпуклая целевая функция, то я бы его для данной задачи не рекомендовал.

  3. 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 есть возможность скачать бесплатную академическую лицензию. Так как я не являюсь участником университетской деятельности, мне этот путь заказан, но согласитесь, подход компаний подкупает, пока не посмотришь на коммерческие ценники.


Перейдём к самому алгоритму.

  1. Шаг

Возьмём симметричную матрицу смежности для некого графа на 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

Нам остаётся только отобразить на графе данные рёбра.

Мы нашли Гамильтонов цикл минимальной длины.

  1. Шаг

Возьмём матрицу крупнее и проведём те же действия

[[ 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

Путь стал немного длиннее чем ранее, но зато решение верно.

  1. Шаг

Для графов большей размерности мы можем проворачивать подобный трюк в цикле до тех пор, пока решение не будет найдено.

Однако даже для такого не серьёзного графа, как этот, 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 дополнительных ограничений, к тому же расчёт затянулся на несколько минут.

Так как с ростом размера входной матрицы количество множеств, на которое может распасться решение возрастает экспоненциально. Нам нужна другая стратегия!

  1. Шаг

Утверждаю, что аналогичного решения можно добиться всего за 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 на предыдущих этапах, должно существовать по крайней мере две связи, соединяющие его с другими подмножествами соответствующего этапа.

Процесс поиска повторять до тех пор, пока не останется только одно множество, оно и будет искомым.


По моим выкладкам получается, что даже в самом худшем случае алгоритму на симметричном полносвязном графе необходимо не больше чем шагов:

\left\lceil \frac{n}{2} \right\rceil+1

и не больше чем дополнительных ограничений:

\sum_{i=6;n\ge 6}^{n}\left\lfloor\sqrt{n}-1\right\rfloor

Число 6 в сумме не случайно, на графах с числом вершин пять и менее дополнительные ограничения совсем не нужны, так как в таком случае не могут появится подциклы. (Для образования подцикла при текущей формулировке задачи требуется минимум три вершины.)

Так для n = 100, нужно не более 523, дополнительных ограничений в худшем случае, а для n=1000, 19613, что вполне не плохо.

Зависимость количества дополнительных ограничений от размера входной матрицы в худшем случае.
Зависимость количества дополнительных ограничений от размера входной матрицы в худшем случае.

Исходя из того, что исходные условие задаются симметричной матрицей смежности, мы можем для экономии памяти отбросить верхнюю или нижнюю её часть, а также главную диагональ.

Уменьшим, таким образом, число переменных для солвера до:

\frac{(n^{2}-n)}{2}

Почему же так важно уменьшить число переменных и дополнительных ограничений? Если под капотом у решателя находится симплекс метод, или его аналог, то формируется матрица:

N*M; N>M

где 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)


  1. Krasnoarmeec
    21.01.2023 11:32
    +4

    Статья супер!
    Хабр опять торт!
    Особенно впечатляет визуализация решений при разном количестве ограничений.
    Сильно жаль, что могу поставить за статью только два плюсика. ????


  1. KongEnGe
    21.01.2023 11:45
    +3

    Очень не хватает формального описания подхода для понимания сути происходящего.


    1. rebuilder Автор
      21.01.2023 12:24
      +1

      Очень не хватает формального описания подхода для понимания сути происходящего.

      Теоретическое обоснование разжевано десятилетия назад. Читать не математику это будет не интересно, а тот к то в теме поймёт опираясь на свой опыт работы с TSP.

      Всё что я сделал, это применил новый подход к выбору ограничивающих условий.


      1. KongEnGe
        21.01.2023 12:46
        +2

        Тем более, если теоретически обосновано, почему не дать ссылку на готовое? А то я к TSP небезразличен, но вот реконструировать подход по косвенным признакам уже лень :)

        То, что я вижу, никак не говорит о том, что в варианте задачи дискретной оптимизации есть уход от NP-полноты. Главное матрицу расстояний попротивнее подогнать на обсчет.


        1. rebuilder Автор
          21.01.2023 13:02

          Понимаете коллега,  где вы тут увидели дискретную оптимизацию? В ней я разочаровался еще, когда делал ветви и границы. Как раз линейное целочисленное программирование, основанное на симплекс методе и алгоритме Гомори, как раз и даёт эти чудесные результаты.


          1. KongEnGe
            21.01.2023 13:08
            +1

            У Вас после жирной "Сути метода" пошла формула как раз из раздела "Дискретная оптимизация":
            https://ru.wikipedia.org/wiki/Задача_коммивояжёра#Формулировка_в_виде_задачи_дискретной_оптимизации

            Поэтому и хотелось увидеть, где случился ранее незамеченный поворот в прорывную методику.


            1. rebuilder Автор
              21.01.2023 13:26

              О! Я так рад, что вы заинтересовались и начали копать. Вы абсолютно правы, за основу была выбрана именно эта формулировка. Но если посмотреть работы Данцига, Фалкерсона и Джонсона, а так же Миллера, Такера и Землина (ссылки на эти работы есть там же в википедии), то можно прийти к следующему выводу:

              Эти авторы пытались засунуть все ограничения в систему уравнений сразу, а так их получается очень много и решатель захлёбывается. Я же предлагаю добавлять только те ограничения, которые нужны по факту в конкретной задаче.

              Надеюсь так яснее. Я думал код приведённый выше сам за себя скажет.


              1. rebuilder Автор
                21.01.2023 15:12

                Тут я не совсем корректно выразился. Выше названные исследователи, конечно же, перебирали комбинации. Если выразить эти комбинации в виде линейных неравенств и отдать их целочисленному решателю, то будет слишком много уравнений.


                1. KongEnGe
                  22.01.2023 11:45

                  Мне просто все происходящее напомнило одну давнюю дискуссию про "коммивояжера за полином" (https://forum.ixbt.com/topic.cgi?id=40:3025)

                  Там была и ссылка на наборы "нехороших" данных для задачи, которая, увы, протухла, но, вроде бы, я нашел именное ее файл в своей домашней помойке. Не хотите ли скормить в свой алгоритм?

                  https://drive.google.com/file/d/11PfpGmwgUCbyg9xeyEeQYpa-niWp53sL


                  1. rebuilder Автор
                    22.01.2023 12:22

                    Не хотите ли скормить в свой алгоритм?

                    Только если вы расскажете, как это прочитать. В архиве что угодно, но только не симметричная матрица смежности.

                    По поводу дискуссии на форуме, как мне показалось, там речь шла про точность самой вводимой матрицы.

                    В тексте статьи был обозначен худший случай, требующий (n/2+1) вызовов решателя. В своих изысканиях специально искал примеры, которые являются адом для решателя, но и там он справлялся не плохо. Хоть и нелегко, но вдвое бил динамическое программирование по размеру матрицы. Если у динамического программирования сложность O(n^2*2^n), то тут я предполагаю она упирается в худшем случае в O(n*2^n) но это прикидочная величина.


                    1. KongEnGe
                      22.01.2023 12:33

                      Не скажу за все, но некоторые наборы симметричные при беглом взгляде (br17). Матрицу вычитывать последовательно, значение за значением, в количестве, определяемом в DIMENSION


              1. guyfawkes
                23.01.2023 23:24

                Простите, но ваш код невозможно читать. Вложенные циклы, переменные с ничего не говорящими названиями. Что происходит, без отладчика совершенно неясно, особенно если человек не понимает вашего хода мысли. Общие моменты текста статьи ещё кое-как маппятся на вашу реализацию, но детали остаются неочевидны.


                1. rebuilder Автор
                  24.01.2023 09:17

                  Признаю код не самый простой, есть грешок. Самые сложные места постарался прокомментировать по ходу. Вложенные циклы можно спокойно развернуть и добавить точки отладчика чтобы посмотреть сложные места. Непонятные переменные только для циклов, для краткости, не люблю, знаете, длинные строки.

                  Если есть конкретные вопросы спрашивайте отвечу. Так долго варился в этом коде что не могу критически оценить его сложность.


            1. thevlad
              21.01.2023 23:36
              +2

              Дискретная оптимизация подразумевает различные методы. В данном случаи используется формулировка, как задачи целочисленного программирования. Прелесть в том, что матрицу подогнать конечно теоретически можно, но на практике это довольно сложно. Уже давал в предыдущей статье автора ссылку на SOTA решатель, основанный на тех же принципах https://en.wikipedia.org/wiki/Concorde_TSP_Solver, он может решать инстансы до десятков тысяч вершин. В то время как другие наивные подходы уже взрываются на паре десятков.


              1. rebuilder Автор
                21.01.2023 23:45

                Благодарю вас за наводку на солвер Concorde, он лучший. Однако после ознакомления с их описание задачи у меня сложилось чёткое убеждение, что они не используют именно целочисленное линейное программирование, скорее там обычное линейное программирование. А так же я не увидел именно точного решения. Возможно, вы мне поможете с этим разобраться?


                1. 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


      1. tessob
        23.01.2023 17:18
        +1

        Было бы интересно взглянуть на то как формулируются линейные ограничения для задачи коммивояжера. Просто из статьи это никак не следует, а сам по себе Симплекс - ни разу не магия.


  1. abutorin
    21.01.2023 12:09

    Из открытых библиотек есть еще GLPK, OR-Tools от Google.


    1. rebuilder Автор
      21.01.2023 12:17
      +1

      GLPK у меня не захотел решать целочисленную задачу, пришлось забраковать.

      А какой решатель встроен OR-Tools, по умолчанию?


      1. abutorin
        21.01.2023 12:24

        А какой решатель встроен OR-Tools, по умолчанию?

        По умолчанию у них свой решатель. Но можно и другие подключить. Правда другие я не пробовал. На линейних задачах встроенный быстрее чем GLPK. В свое время проводил сравнение SciPy. GLPK, OR-Tools по скорости и максимальной размерности задачи. SciPy отказался решать систему с 5000 переменных. GLPK и OR-Tools с ними справляются, последний несколько быстрее.


        1. rebuilder Автор
          21.01.2023 12:36

          КДПВ для статьи сформирована SciPy и для неё потребовалось 280875 переменных, нормально справился. Обратите внимание, что в последних версиях scipy.optimize.linprog завезли классные быстрые решатели. Особенно меня порадовал метод глубинной точки, для нецелочисленных задач.


          1. abutorin
            21.01.2023 13:04

            для неё потребовалось 280875 переменных

            Тут наверное пора уточнить что такое "переменная". Я подрузумелвал задачу оптимизации матрицы в которой больше 5000 строк и столбцов, т.е. там 25 млн. ячеек.

            Но спасибо за наводку, посмотрю свежую версию SciPy.


            1. rebuilder Автор
              22.01.2023 15:49

              Стандартный решатель  OR-Tools, меня совсем не порадовал, по скорости он на последнем месте.


  1. DonAgosto
    21.01.2023 15:31

    На печатной плате нужно просверлить монтажные отверстия (произвести пайку / прочее действие)
    Заметная польза будет, если нахождение в «вершине» дешевле, чем перемещение по «ребру». К сожалению чаще бывает наоборот.


    1. rebuilder Автор
      21.01.2023 16:55

      Всё так, но ведь и операции могут быть разными. Например, монтажный пистолет работает достаточно шустро, прозвонка дорожек мультиметром то же весьма скоростная  операция. К сожалению, не владею реальными кейсами.


  1. NiceITMan
    21.01.2023 22:39

    Отличная статья да и люди тут важные )


    1. rebuilder Автор
      21.01.2023 22:39

      Да не, мы тут все люди простые, вот автор статьи обычный разраб. Просто у него хобби искать решения не решаемых задач.


  1. tzlom
    22.01.2023 12:07

    Примерно так работает Active Set Solver для QP проблем.

    Правда в нем еще есть шаг релаксации ограничений, т.к. можно в процессе построить ограничения которые будут загонять в локальный минимум.

    Кстати для AS солвера мне помогло не решать проблему полностью заново, а дорешивать её для новых условий , сложность падает с n^3 до ~n^2

    Интересно, что будет если уточняющие ограничения вводить с условием >=1 ?


    1. rebuilder Автор
      22.01.2023 12:37

      Интересно, что будет если уточняющие ограничения вводить с условием >=1 ?

      Решение всё равно найдётся но, большим числом неравенств


      1. tzlom
        22.01.2023 12:41

        Это теория или попробовали? Ограничение на замкнутость должно дополнять это ограничение в любом случае, разве нет?


        1. rebuilder Автор
          22.01.2023 12:47
          +1

          Попробовал, множества всё равно связываются хоть и одной связью. За счёт повторов связывания, рано или поздно стягивает всю конструкцию в монолит. Но получается медленнее чем с двумя связями.