С 1 по 7 ноября проходила квалификация Yandex Cup 2022. В секции Алгоритм: Марафон организаторы предложили интересную задачу программирования в ограничениях, обобщения известной задачи коммивояжёра, задачу поиска маршрута (vehicle routing problem). В статье расскажу о своем решении на основе Google OR-tools.

Пример маршрутов для первой задачи

Задача состояла в нахождении маршрутов перевоза самокатов для нескольких автомобилей от места, где самокат оставили, до парковки. Автомобиль имеет ограниченный запас хода и ограниченную вместимость. Учитывая сложность подобного рода задач и ограничение по времени (неделя по условиям конкурса), сразу начал смотреть в сторону доступных солверов.

Первым в поиске попался pycsp. Подробно останавливаться не буду. Отмечу понятный интерфейс и удобство описания задачи. Решение было написано и протестировано на модельной задаче за пару часов. С реальными задачами из конкурса, где 1000 узлов против 10 в модельной, солвер не справился. Трансляция в xml не закончилась за ночь. Сделал xml (~5Мб) с помощью jinja, запустил ACE напрямую. Солвер съел всю память(~12Гб), но не смог её прочесть в себя. На этом я с ним попрощался. Потратив первые два дня из выделенных.

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

Составлял решение последовательно модифицируя один из примеров. Всё начинается с описания расстояний между узлами. 

from ortools.constraint_solver import pywrapcp as cp
# nodes - количество узлов
# cars - количество машин
# startn_node = 0 для всех задач - узел отправления автомобилей
# D[i, j] - матрица расстояний из i в j
manager = cp.RoutingIndexManager(nodes, cars, start_node)
routing = cp.RoutingModel(manager)

def distance_callback(from_index, to_index):
    from_node = manager.IndexToNode(from_index)
    to_node = manager.IndexToNode(to_index)
    return D[from_node, to_node]

distance_callback_index = routing.RegisterTransitCallback(distance_callback)
# Define cost of each arc.
routing.SetArcCostEvaluatorOfAllVehicles(distance_callback_index)

Дальше необходимо задать ограничение на запас хода. Так же берем из примера.

routing.AddDimensionWithVehicleCapacity(
    distance_callback_index,
    0,  # no slack
    car_limits,  # vehicle maximum travel distance
    True,  # start cumul to zero
    'Distance'
)

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

def capacity_callback(from_index):
    from_node = self._manager.IndexToNode(from_index)
    return bike_cost[from_node]  # 1 - самокат, -1 - парковка

capacity_callback_index = routing.RegisterUnaryTransitCallback(capacity_callback)

routing.AddDimensionWithVehicleCapacity(
    capacity_callback_index,
    0,  # no slack
    [capacity] * cars,
    True,  # start cumul to zero
    'Capacity'
)

d = routing.GetDimensionOrDie('Capacity')
for car in range(cars):
    end_index = routing.End(car)
    solver.Add(d.CumulVar(end_index) == 0)

Очевидно, все узлы посетить невозможно. Даём возможность пропустить любой узел, но со штрафом. Размер штрафа в целевую функцию выбрал заведомо больше, чем длина маршрута, COEF_NODE = 1 000 000.

for node in range(1, nodes):
    routing.AddDisjunction([manager.NodeToIndex(node)], COEF_NODE)

По условию задачи автомобиль не обязан возвращаться в стартовый узел, но OR-tools так не умеет. Машины всегда будут возвращаться в начало. В документации есть подсказка, как добиться нужного поведения: нужно добавить доп узлы с нулевым расстоянием из любого другого. Расставляем 0 и inf в матрицу так, чтоб автомобиль всегда проходил через промежуточный узел и только один.

n = D.shape[0]
D_new = np.zeros((n + cars, n + cars), dtype=D.dtype)
D_new[:n, :n] = D  # основная матрица без изменений
D_new[:n, 0] = inf  # возвращения только через промежуточный узел
D_new[n:, n:] = inf  # между этими узлами нельзя передвигаться
D_new[-cars:, 1:n] = inf  # из промежуточного узла только в 0
D_new -= np.diag(np.diag(D_new, 0))

bike_cost = bike_cost + [0, ] * cars

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

search_parameters = cp.DefaultRoutingSearchParameters()
search_parameters.use_full_propagation = True
solution = routing.SolveWithParameters(search_parameters)

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

Попробовал ослабить ограничения, вместо жёстких (hard) сделать мягкие (soft), когда за превышение идет штраф, но ограничения как такового нет. Получилось не с первой попытки, но завелось. Размер штрафа пытался подобрать исходя из лучших результатов первой задачи.

# запас хода
routing.AddDimensionWithVehicleCapacity(
    distance_callback_index,
    0,  # no slack
    [100 * lim for lim in car_limits],  # vehicle maximum travel distance
    True,  # start cumul to zero
    'Distance'
)

d = routing.GetDimensionOrDie('Distance')
for car, limit in enumerate(car_limits):
    index = self._routing.End(car)
    d.SetCumulVarSoftUpperBound(index, limit, COEF_DISTANCE)

# вместимость
routing.AddDimensionWithVehicleCapacity(
    capacity_callback_index,
    0,  # no slack
    [nodes] * cars,
    False,  # start cumul to zero
    'Capacity'
)
d = routing.GetDimensionOrDie('Capacity')
for n in range(nodes):
    index = manager.NodeToIndex(n)
    d.CumulVar(index).SetMin(-nodes)
    d.SetCumulVarSoftLowerBound(index, 0, COEF_CAPACITY)
    d.SetCumulVarSoftUpperBound(index, capacity, COEF_CAPACITY)

for car in range(cars):
    end_index = routing.End(car)
    d.SetCumulVarSoftUpperBound(end_index, 0, COEF_NODE)  # they will be discarded

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

  1. ищем маршруты без ограничения на вместимость

  2. улучшаем маршруты с мягким ограничением на вместимость

  3. улучшаем маршруты с жёстким ограничением на вместимость

Но даже такой подход не прошел порог в 10 000 самокатов в сумме из 30 задач, остановившись на 9 300. Полное решение здесь.

Заключение

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