С вами Алексей Ложкинс, эксперт по анализу данных и машинному обучению в ПГК Диджитал. Мы разрабатываем цифровые продукты для логистической отрасли, в первую очередь, для ж/д перевозок.
В кулуарах московского офиса ПГК мы обсуждаем и нерабочие темы. Топовую строчку в темах неформального общения занимает отпуск. Мы решили рассмотреть задачу планирования отпуска, как задачу оптимизации маршрута по выбранным достопримечательностям. Для этого воспользовались классической постановкой задачи коммивояжера.
Моделирование маршрута в виде задачи коммивояжера позволит построить маршрут по всем запланированным локациям без повторений с заданным критерием качества (время, стоимость). Рассмотрим несколько подходов к решению оптимизационной задачи (TSP) с использованием пакета ORTools.
О задаче коммивояжера (TSP)
Доставка разных товаров очень популярная услуга сегодня. Чтобы снизить себестоимость и цену на услуги для потребителя, заказы объединяют и доставляются одним курьером. Здесь появляется так называемый эффект масштаба, когда себестоимость доставки единицы заказа снижается с увеличением их количества.
При рассмотрении одного курьера и пула заказов, которые он должен доставить, возникает классическая задача коммивояжера. В более сложном случае, когда у нас несколько курьеров, задача коммивояжера эволюционирует в vehicle routing problem, но это уже другая история.
Здесь рассмотрим решение задачи коммивояжера для построения туристического маршрута на примере локаций на территории России. Требования к маршруту могут быть разными, рассмотрим наиболее популярные: самый быстрый или самый дешевый.
Туристические локации: данные
Куда мы хотим поехать? Соберем пожелания в список. Дарья Васильева (координатор управления анализа данных и машинного обучения)подготовила свой вариант такого списка, с которым можно ознакомиться здесь, отдельное спасибо ей за это! Помимо локаций потребуются ребра - связь между локациями.
Рассмотрим неполный граф (возможны не все соединения между городами, а наиболее рациональные ребра). География предлагаемых мест достаточно обширная, и время следования ж/д или авто транспортом может занимать рабочую неделю, поэтому сделали акцент на авиаперевозки. Если хотите насытить поездку ж/д антуражем, то достаточно наполнить входные данные тарифами и временем в ходу для интересующих ребер. Модель можно использовать и для грузоперевозок по железной дороге в том числе. С графом допустимых перемещений можно ознакомиться здесь.
Технический стэк
Прежде чем переходить к технической части, обсудим инструменты и подход. Математический аппарат для моделирования задачи - целочисленное линейное программирование (ILP), относится к точным методам. Постановка в ILP позволяет применить уже готовые пакеты (солверы) для решения оптимизационной задачи и получить точное решение на выходе.
Солверы можно разделить на бесплатные open-source солверы (cbc, scip, highs, glpk и др.) и коммерческие (gurobi, copt, cplex, optverse и др.). Так как солверы применяются только для моделирования и решения задач, то спектр отличий у них не большой: API и производительность. Вопрос API решается сторонними обертками над солверами, которые позволяют моделировать и подключать различные солверы для решения задач (например, PuLP, Pyomo, ORTools, FeLooPy, GAMS и др.). Как следствие, коммерческие солверы "почти везде" превосходят по производительности open-source, больше фактов здесь.
Зафиксируем стэк. Реализацию модели сделаем на python, среда для моделирования ORTools, солвер open-source cp-sat (ORTools), визуализацию будем проводить в plotly, немного поиспользуем pandas и networkx.
По поводу солвера, cp-sat не совсем решатель задач целочисленного линейного программирования, но в том числе покрывает этот класс задач. Детально в специфику не будем уходить (ключевые слова для тех кто хочет: constraint programming, boolean satisfability). Кроме того, ORTools имеет отдельный набор алгоритмов для решения задачи vehicle routing problem, которая в частном случае эквивалентна задаче TSP. Воспользуемся этими алгоритмами для валидации результата модели ILP.
Моделирование
Воспользуемся симметричностью матрицы переходов между узлами и сформулируем минималистичную модель ILP задачи планирования маршрута на условиях TSP. Введем обозначения, ограничения и целевую функцию:
- множество узлов сети;
- расстояние/затраты/время движения по ребру ;
- бинарная переменная, индикатор выбора ребра в маршрут (=1 если ребро в маршруте, 0 - в противном случае);
Начнем разбор модели с целевой функции - традиционно это минимизация длины маршрута, но как вариации могут рассматриваться минимизация затрат или минимизация времени в пути, другие варианты тоже возможны. Ограничение (1) обязывает модель организовать один вход и один выход. Для каждого подмножества создаем ограничения типа (2), они обеспечивают создание связного маршрута. Ограничение симметрии (3), неважно, в каком направлении идем по ребру, важно что идем (матрица переходов симметричная). Ограничения (4) - бинарность переменных.
Комментарии к ограничениям. Ограничения типа (3) сократятся на presolve этапе солвера. Можно сформулировать модель без этого ограничения, но запись модели будет немного сложнее. Количество ограничений типа (2) растет экспоненциально в зависимости от , что выглядит не очень практично. С этим поборемся далее.
Python реализация модели
Построим модель без ограничения (2), посмотрим, что означает связность маршрута на практике. В качестве целевой функции будем использовать минимизацию затрат на трансфер между локациями.
Реализация упрощенной версии модели без ограничений типа (2)
# Упрощенная версия модели без ограничений типа (2)
from ortools.sat.python import cp_model
import pandas as pd
import numpy as np
# Считаем входные данные
df_dist = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_dist.csv", sep=";", encoding="cp1251")
df_location = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_coord.csv", sep=";", encoding="cp1251")
obj_type = "cost" # Тип целевой функции. Допустимые значения ["cost", "travel_time"]
m = cp_model.CpModel()
# Инициализация переменных
df_dist["var_id"] = np.ogrid[:df_dist.shape[0]]
df_dist["var"] = df_dist["var_id"].apply(lambda x: m.NewBoolVar(f"x_{x}"))
# Создаем ограничения типа (1)
dct_constr_1 = df_dist.groupby("to")["var"].sum().to_dict()
for from_node, var_sum in dct_constr_1.items():
m.Add(var_sum == 2)
# Создаем ограничения типа (3)
dct_constr_3 = df_dist.set_index(["from", "to"]).to_dict()["var"]
for (from_node, to_node), var_ in dct_constr_3.items():
m.Add(var_ == dct_constr_3[to_node, from_node])
m.Minimize(sum(df_dist[obj_type] * df_dist["var"]))
# Инициализация solver
solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
# Решение задачи
status = solver.Solve(m)
# Проверяем статус
if status == cp_model.OPTIMAL:
print("-" * 40)
gap_val = abs((solver.ObjectiveValue() - solver.BestObjectiveBound()) / max(1, solver.ObjectiveValue()))
# Извлекаем результат
df_dist["sol"] = df_dist["var"].apply(lambda x: solver.Value(x))
df_dist_opt = df_dist[df_dist["sol"] > 0]
print(f"-- Found solution {solver.ObjectiveValue()}")
print(f"-- GAP: {round(gap_val, 2)}%")
print(f"-- Solved in {solver.WallTime()} sec.")
print(f"-- Costs {df_dist_opt['cost'].sum() / 2} rub")
print(f"-- Duration {df_dist_opt['travel_time'].sum() / 2} h")
Функция визуализация решения в plotly
# Нанесение результата на карту
import plotly.graph_objects as go
def plot_route(df_graph: pd.DataFrame):
# Начало и окончание маршрута в МСК
df_location["color"] = np.where(df_location["id"] == 0, "red", "blue")
df_location["size"] = np.where(df_location["id"] == 0, 14, 10)
df_location['text'] = df_location["id"].astype(str) + ", " + df_location["name"]
# Добавляем координаты и наименование точек отправления и прибытия
df_loc_from = df_location.copy()
df_loc_from.columns = [col + "_from" for col in df_loc_from.columns]
df_loc_to = df_location.copy()
df_loc_to.columns = [col + "_to" for col in df_loc_to.columns]
df_travel_map = df_graph.merge(df_loc_from, how="left", left_on="from",
right_on="id_from")
df_travel_map = df_travel_map.merge(df_loc_to, how="left", left_on="to",
right_on="id_to")
# Наносим вершины
fig = go.Figure(data=go.Scattergeo(
lon = df_location['lon'],
lat = df_location['lat'],
text = df_location['text'],
mode = 'markers',
name = "Локации",
marker = dict(
size = df_location["size"],
opacity = 0.8,
color = df_location["color"],
),
))
# Наносим ребра
for row in df_travel_map.itertuples():
fig.add_trace(
go.Scattergeo(
lon = [row.lon_from, row.lon_to],
lat = [row.lat_from, row.lat_to],
mode = 'lines',
line = dict(width = 1,color = 'blue'),
opacity = 0.8,
name=row.name_from + " - " + row.name_to
),
)
fig.update_layout(
title = 'Travel route in Russia <br>',
geo = dict(
showland = True,
landcolor = "rgb(250, 250, 250)",
subunitcolor = "rgb(217, 217, 217)",
countrycolor = "rgb(217, 217, 217)",
countrywidth = 0.5,
subunitwidth = 0.5
),
)
fig.show()
С целевой функцией минимизации затрат получили оптимальное решение задачи с общими затратами 93 876 руб (115.05 часов). Параметр obj_type
позволяет переключить целевую функцию на минимизацию времени в пути, тогда оптимальное решение будет 70.9 часов (131 493 руб). Нанесем результат на карту для целевой функции минимизации затрат.
Не такой результат мы с вами ожидали, маршрут получился разбитым на части - несвязным. Ограничения типа (2) устраняют эту проблему, но цена достаточно высока - громаднейшая модель. Предложение - решать проблемы по мере их поступления.
Исходя из решения задачи, мы можем определить несвязные подграфы и специально для них создать ограничения типа (2) - добавим заплатку. Само ограничение "разрывает" локальные циклы, обязывая создавать хотя бы одну связь с внешним множеством вершин. Новые ограничения приводят к изменению модели, следовательно, необходимо искать новое решение.
В некоторых солверах есть возможность создавать ленивые ограничения (lazy constraints). Это ограничения, которые изначально находятся в неактивном состоянии, но при их нарушении или по другим критериям активации они добавляются в модель непосредственно в процессе решения задачи.
Возвращаясь к нашей задаче. Было бы неплохо лениво создавать ограничения типа (2) при получении решения с нарушением связности. Но, к сожалению, ortools не имеет примеров или хорошей документации метода set_is_lazy, поэтому мы воспользуемся циклами python. Связные подграфы будем определять с помощью библиотеки networkx. Критерий остановки - получен связный маршрут.
*Работу ленивых ограничений можно посмотреть у Gurobi, не реклама.
Алгоритм итеративного добавления ограничения типа (2)
# Итеративное добавление ограничений типа (2)
from ortools.sat.python import cp_model
import networkx as nx
import pandas as pd
import numpy as np
# Считаем входные данные
df_dist = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_dist.csv", sep=";", encoding="cp1251")
df_location = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_coord.csv", sep=";", encoding="cp1251")
obj_type = "travel_time" # Тип целевой функции. Допустимые значения ["cost", "travel_time"]
m = cp_model.CpModel()
# Инициализация переменных
df_dist["var_id"] = np.ogrid[:df_dist.shape[0]]
df_dist["var"] = df_dist["var_id"].apply(lambda x: m.NewBoolVar(f"x_{x}"))
# Создаем ограничения типа (1)
dct_constr_1 = df_dist.groupby("to")["var"].sum().to_dict()
for from_node, var_sum in dct_constr_1.items():
m.Add(var_sum == 2)
# Создаем ограничения типа (3)
dct_constr_3 = df_dist.set_index(["from", "to"]).to_dict()["var"]
for (from_node, to_node), var_ in dct_constr_3.items():
m.Add(var_ == dct_constr_3[to_node, from_node])
m.Minimize(sum(df_dist[obj_type] * df_dist["var"]))
# Инициализация solver
solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
# Решение задачи
status = solver.Solve(m)
# Извлекаем результат
df_dist["sol"] = df_dist["var"].apply(lambda x: solver.Value(x))
df_dist_opt = df_dist[df_dist["sol"] > 0]
# ----------- Итерации по добавлению ограничений типа (2). НАЧАЛО ------------------
# Кол-во связных подграфов
route_tmp = nx.from_pandas_edgelist(df_dist_opt, "from", "to")
lst_sg = list(nx.connected_components(route_tmp)) # Список узлов в связных подграфах
# Для статистики
cnt_constrs = 0 # Кол-во добавленных ограничений типа (2)
cnt_iter = 0 # Кол-во итераций решения задачи ILP
solver_time = 0 # Суммарное время решения задач солвером
# Повторять, пока имеем более одного несвязного подграфа
while len(lst_sg) > 1:
cnt_iter += 1
for sg in lst_sg:
cnt_constrs += 1
# Отбираем необходимые переменные для суммирования
filt_i = df_dist["from"].isin(sg)
filt_j = df_dist["to"].isin(sg)
filt_i_less_j = df_dist["from"] < df_dist["to"]
filt_all = filt_i_less_j & filt_i & filt_j
var_sum = df_dist[filt_all]["var"].sum()
# Добавляем ограничение типа (2) в модель
m.Add(var_sum <= len(sg) - 1)
status = solver.Solve(m)
solver_time += solver.WallTime()
# Извлекаем результат
df_dist["sol"] = df_dist["var"].apply(lambda x: solver.Value(x))
df_dist_opt = df_dist[df_dist["sol"] > 0]
# Кол-во связных подграфов
route_tmp = nx.from_pandas_edgelist(df_dist_opt, "from", "to")
lst_sg = list(nx.connected_components(route_tmp)) # Список узлов в связных подграфах
# ----------- Итерации по добавлению ограничений типа (2). КОНЕЦ ------------------
# Проверяем статус
if status == cp_model.OPTIMAL:
print("-" * 40)
gap_val = abs((solver.ObjectiveValue() - solver.BestObjectiveBound()) / max(1, solver.ObjectiveValue()))
print(f"-- Found solution {solver.ObjectiveValue()}")
print(f"-- GAP: {round(gap_val, 2)}%")
print(f"-- Solved in {solver_time} sec.")
print(f"-- Constraints (2) added {cnt_constrs}")
print(f"-- Iterations count {cnt_iter}")
print(f"-- Costs {df_dist_opt['cost'].sum() / 2} rub")
print(f"-- Duration {df_dist_opt['travel_time'].sum() / 2} h")
Полная версия исходной модели (все ограничения типа (2))
# Полная модель - инициализируем все ограничения типа (2)
from ortools.sat.python import cp_model
import pandas as pd
import numpy as np
from itertools import chain, combinations
# Считаем входные данные
df_dist = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_dist.csv", sep=";", encoding="cp1251")
df_location = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_coord.csv", sep=";", encoding="cp1251")
obj_type = "travel_time" # Тип целевой функции. Допустимые значения ["cost", "travel_time"]
m = cp_model.CpModel()
# Инициализация переменных
df_dist["var_id"] = np.ogrid[:df_dist.shape[0]]
df_dist["var"] = df_dist["var_id"].apply(lambda x: m.NewBoolVar(f"x_{x}"))
# Создаем ограничения типа (1)
dct_constr_1 = df_dist.groupby("to")["var"].sum().to_dict()
for from_node, var_sum in dct_constr_1.items():
m.Add(var_sum == 2)
# Создаем ограничения типа (2)
def powerset(iterable):
s = list(iterable)
return chain.from_iterable(combinations(s, r) for r in range(len(s)+1) if 1 < r < len(s) - 1)
cnt_constrs = 0
for lst_nodes in powerset(df_location["id"]):
cnt_constrs += 1
# Отбираем необходимые переменные для суммирования
filt_i = df_dist["from"].isin(lst_nodes)
filt_j = df_dist["to"].isin(lst_nodes)
filt_i_less_j = df_dist["from"] < df_dist["to"]
filt_all = filt_i_less_j & filt_i & filt_j
var_sum = df_dist[filt_all]["var"]
if len(var_sum) != 0:
cnt_constrs += 1
var_sum = var_sum.sum()
# Добавляем ограничение типа (2) в модель
m.Add(var_sum <= len(lst_nodes) - 1)
# Создаем ограничения типа (3)
dct_constr_3 = df_dist.set_index(["from", "to"]).to_dict()["var"]
for (from_node, to_node), var_ in dct_constr_3.items():
m.Add(var_ == dct_constr_3[to_node, from_node])
m.Minimize(sum(df_dist[obj_type] * df_dist["var"]))
# Инициализация solver
solver = cp_model.CpSolver()
solver.parameters.log_search_progress = True
# Решение задачи
status = solver.Solve(m)
# Проверяем статус
if status == cp_model.OPTIMAL:
print("-" * 40)
gap_val = abs((solver.ObjectiveValue() - solver.BestObjectiveBound()) / max(1, solver.ObjectiveValue()))
print(f"-- Found solution {solver.ObjectiveValue()}")
print(f"-- GAP: {round(gap_val, 2)}%")
print(f"-- Solved in {solver.WallTime()} sec.")
print(f"-- Constraints (2) added {cnt_constrs}")
print(f"-- Costs {df_dist_opt['cost'].sum() / 2} rub")
print(f"-- Duration {df_dist_opt['travel_time'].sum() / 2} h")
# Извлекаем результат
df_dist["sol"] = df_dist["var"].apply(lambda x: solver.Value(x))
df_dist_opt = df_dist[df_dist["sol"] > 0]
Что изменилось в результате? Полученный маршрут представляет собой один большой цикл без повторения вершин - гамильтонов цикл. Это результат, к которому стремились! Берем на заметку и отправляемся в путь!
Проверим решение, для этого воспользуемся готовым модулем от ORTools для решения задачи VRP. Потребуется подстроить входные данные под алгоритм. Матрица переходов между узлами у нас неполная, поэтому воспользуемся сравнительно большими числами в качестве заглушек для пропущенных значений. Количество единиц транспорта = 1 в случае задачи TSP, узел начала и окончания маршрута - Москва (0). Отмечу, что матрица "расстояний" должна быть целочисленной.
Пакетная версия решения задачи VRP от ORTools
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
# Считаем входные данные
df_dist = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_dist.csv", sep=";", encoding="cp1251")
df_location = df_cab = pd.read_csv("https://raw.githubusercontent.com/Lozkins/mos/master/data/06_tsp_coord.csv", sep=";", encoding="cp1251")
obj_type = "cost" # Тип целевой функции. Допустимые значения ["cost", "travel_time"]
large_M = 10**5
round_mult = 1
def create_data_model():
"""
Подготовка данных к запуску
"""
data = {}
dct_dist = df_dist.set_index(["from", "to"]).to_dict()[obj_type]
# Матрица расстояний
data["distance_matrix"] = [
[int(dct_dist.get((i, j), 0 if i == j else large_M) * round_mult) for j in range(df_location.shape[0])]
for i in range(df_location.shape[0])]
# Кол-во транспортных средств - для задачи TSP - 1
data["num_vehicles"] = 1
# Вершина начала и завершения маршрута
data["depot"] = 0
return data
def print_solution(manager, routing, solution):
"""
Вывод решения
"""
print(f"Значение целевой функции: {solution.ObjectiveValue()} {obj_type}")
index = routing.Start(0)
plan_output = "Маршрут:\n"
route_distance = 0
while not routing.IsEnd(index):
plan_output += f" {manager.IndexToNode(index)} ->"
previous_index = index
index = solution.Value(routing.NextVar(index))
route_distance += routing.GetArcCostForVehicle(previous_index, index, 0)
plan_output += f" {manager.IndexToNode(index)}\n"
plan_output += f"Параметр маршрута: {route_distance} {obj_type}\n"
def get_routes(solution, routing, manager):
"""
Извлечение маршрута из модели
"""
routes = []
for route_nbr in range(routing.vehicles()):
index = routing.Start(route_nbr)
route = [manager.IndexToNode(index)]
while not routing.IsEnd(index):
index = solution.Value(routing.NextVar(index))
route.append(manager.IndexToNode(index))
routes.append(route)
return routes
def main():
"""
Запуск алгоритма
"""
# Инициализация данных
data = create_data_model()
# Создание менеджера маршрутизации
manager = pywrapcp.RoutingIndexManager(
df_location.shape[0], data["num_vehicles"], data["depot"])
# Инициализация модели маршрутизации
routing = pywrapcp.RoutingModel(manager)
def distance_callback(from_index, to_index):
from_node = manager.IndexToNode(from_index)
to_node = manager.IndexToNode(to_index)
return data["distance_matrix"][from_node][to_node]
transit_callback_index = routing.RegisterTransitCallback(distance_callback)
# Определяем вес ребра (стоимость или время в пути)
routing.SetArcCostEvaluatorOfAllVehicles(transit_callback_index)
# Инициализация начального решений
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
routing_enums_pb2.FirstSolutionStrategy.PATH_CHEAPEST_ARC
)
# Solve the problem.
solution = routing.SolveWithParameters(search_parameters)
# Print solution on console.
if solution:
print_solution(manager, routing, solution)
routes = get_routes(solution, routing, manager)
# Display the routes.
for i, route in enumerate(routes):
print('Route', i, route)
return route
route = main()
# Преобразуем в таблицу
lst_pairs = list(zip(route[:-1], route[1:]))
df_route_ortools = pd.DataFrame(lst_pairs, columns=["from", "to"])
df_route_ortools = df_dist.merge(df_route_ortools, how="right", on=["from", "to"])
Сведем результаты экспериментов в одну таблицу. Представлю три версии модели: полная модель со всеми ограничениями типа (2) (full), итеративное добавление ограничения (2) (iter) и пакетная версия решения задачи VRP от ORTools (ortools). Замеряемые параметры: время расчета, кол-во ограничений типа (2), кол-во итераций решения ILP, затраты на маршруте и длительность движения по маршруту.
Резюме
Теоретическая постановка задачи со всеми ограничениями типа (2) оказалась наименее производительной, с точки зрения времени поиска оптимального решения. Итеративная постановка подразумевает создание ограничений по мере их необходимости, как результат - создание до 10 ограничений типа (2) вместо ~131к ограничений на нашем сценарии. Масштабирование эффекта итераций на другие наборы данных не работает. Остается вероятность сценария, когда будут генерироваться значительно бОльшие объемы ограничений типа (2).
Эксперименты проводили на сравнительно малом кол-ве вершин, поэтому сравнивать итеративную постановку и ortools будет нерепрезентативно. Ожидание - ortools будет быстрее.
Так куда же ехать? У нас была возможность построить самый быстрый и самый дешевый маршрут. Готовые маршруты по туристическим местам от Даши Васильевой:
Самый дешевый маршрут (101 833 руб / 124,3 ч): Москва -> Волгоград -> Карелия/Петрозаводск -> Псков -> Санкт-Петербург -> Село Териберка/Мурманск -> Калининград -> Листвянка/Иркутск -> Красноярск -> Южно-Сахалинск -> Благовещенск -> Шерегеш/Новокузнецк -> Новосибирск -> Алтай/Горно-Алтайск -> Казань -> Сочи -> Москва
Самый быстрый маршрут (137 770 руб / 71,55 ч): Москва -> Село Териберка/Мурманск -> Калининград -> Псков -> Карелия/Петрозаводск -> Волгоград -> Санкт-Петербург -> Казань -> Сочи -> Листвянка/Иркутск -> Красноярск -> Благовещенск -> Южно-Сахалинск -> Алтай/Горно-Алтайск -> Новосибирск -> Шерегеш/Новокузнецк -> Москва
Будем рады за предложения по локациям и критериям. Быть может, кто-то из нас решится на длительное турне и обкатает решение вживую. Техническая версия статьи здесь