Задача коммивояжера (TSP) точное решение — метод целочисленного линейного программирования (Integer programming)

#------------------------------------------------------------------------------------
# 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_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()

#------------------------------------------------------------------------------------
# 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)