Документация Engee

Задача о коммивояжере

Это руководство было изначально предоставлено Дэниелом Шермером (Daniel Schermer).

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

В нем используются следующие пакеты.

using JuMP
import GLPK
import Random
import Plots

Математическая формулировка

Предположим, дан полный граф , где  — это множество вершин (или городов), а  — множество ребер (или дорог). Для каждой пары вершин с ребром связан вес (или расстояние) .

В этом руководстве предполагается, что задача симметрична, то есть .

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

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

Целевая функция

Цель состоит в том, чтобы минимизировать длину маршрута (из-за предполагаемой симметрии вторая сумма содержит только ):

Обратите внимание, что вместо этого можно также использовать следующую целевую функцию:

Ограничения

В нашей формулировке имеется четыре класса ограничений.

Во-первых, ввиду предполагаемой симметрии должны выполняться следующие ограничения:

Во-вторых, для каждой вершины необходимо выбрать ровно два ребра, соединяющих ее с другими вершинами графа :

В-третьих, не допускаются циклы:

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

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

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

Реализация

В JuMP есть два способа устранения подциклов, и оба они будут показаны ниже:

  • итеративное решение новой модели, которая включает ранее выявленные подциклы;

  • добавление подциклов в виде отложенных ограничений.

Данные

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

function generate_distance_matrix(n; random_seed = 1)
    rng = Random.MersenneTwister(random_seed)
    X = 100 * rand(rng, n)
    Y = 100 * rand(rng, n)
    d = [sqrt((X[i] - X[j])^2 + (Y[i] - Y[j])^2) for i in 1:n, j in 1:n]
    return X, Y, d
end

n = 40
X, Y, d = generate_distance_matrix(n)
([23.603334566204694, 34.651701419196044, 31.27069683360675, 0.790928339056074, 48.86128300795012, 21.096820215853597, 95.1916339835734, 99.99046588986135, 25.166218303197184, 98.66663668987997  …  46.33502459235987, 18.582130997265377, 11.198087695816717, 97.6311881619359, 5.161462067432709, 53.80295812064833, 45.56920516275036, 27.93951106725605, 17.824610354168602, 54.89828719625274], [37.097066286146884, 89.41659192657593, 64.80537482231894, 41.70393538841062, 14.456554241360564, 62.24031828206811, 87.23344353741976, 52.49746566167794, 24.159060827129643, 88.48369255734127  …  66.12321555087209, 19.45678064479248, 39.3193497656424, 99.07406554003964, 55.03342139580574, 58.07816346526631, 76.83586278313636, 51.952465724186084, 51.486297110544356, 99.81360570779374], [0.0 53.473350122820904 … 15.506244459460921 70.09092934998034; 53.473350122820904 0.0 … 41.49527995497558 22.760099542720535; … ; 15.506244459460921 41.49527995497558 … 0.0 60.9096566304971; 70.09092934998034 22.760099542720535 … 60.9096566304971 0.0])

Для модели JuMP сначала инициализируем объект модели. Затем создадим двоичные переменные решения и добавим целевую функцию и ограничения. Определяя матрицу x как Symmetric, мы избегаем необходимости явным образом добавлять ограничения, согласно которым x[i, j] == x[j, i].

function build_tsp_model(d, n)
    model = Model(GLPK.Optimizer)
    @variable(model, x[1:n, 1:n], Bin, Symmetric)
    @objective(model, Min, sum(d .* x) / 2)
    @constraint(model, [i in 1:n], sum(x[i, :]) == 2)
    @constraint(model, [i in 1:n], x[i, i] == 0)
    return model
end
build_tsp_model (generic function with 1 method)

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

function subtour(edges::Vector{Tuple{Int,Int}}, n)
    shortest_subtour, unvisited = collect(1:n), Set(collect(1:n))
    while !isempty(unvisited)
        this_cycle, neighbors = Int[], unvisited
        while !isempty(neighbors)
            current = pop!(neighbors)
            push!(this_cycle, current)
            if length(this_cycle) > 1
                pop!(unvisited, current)
            end
            neighbors =
                [j for (i, j) in edges if i == current && j in unvisited]
        end
        if length(this_cycle) < length(shortest_subtour)
            shortest_subtour = this_cycle
        end
    end
    return shortest_subtour
end
subtour (generic function with 1 method)

Объявим вспомогательную функцию selected_edges(), которая будет неоднократно использоваться в дальнейшем.

function selected_edges(x::Matrix{Float64}, n)
    return Tuple{Int,Int}[(i, j) for i in 1:n, j in 1:n if x[i, j] > 0.5]
end
selected_edges (generic function with 1 method)

Другие вспомогательные функции для вычисления подциклов:

subtour(x::Matrix{Float64}) = subtour(selected_edges(x, size(x, 1)), size(x, 1))
subtour(x::AbstractMatrix{VariableRef}) = subtour(value.(x))
subtour (generic function with 3 methods)

Итеративный метод

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

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

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

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

 
iterative_model = build_tsp_model(d, n)
optimize!(iterative_model)
@assert is_solved_and_feasible(iterative_model)
time_iterated = solve_time(iterative_model)
cycle = subtour(iterative_model[:x])
while 1 < length(cycle) < n
    println("Found cycle of length $(length(cycle))")
    S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
    @constraint(
        iterative_model,
        sum(iterative_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
    )
    optimize!(iterative_model)
    @assert is_solved_and_feasible(iterative_model)
    global time_iterated += solve_time(iterative_model)
    global cycle = subtour(iterative_model[:x])
end

objective_value(iterative_model)
525.7039004442727
time_iterated
0.03748679161071777

Для быстрой проверки работоспособности мы визуализируем оптимальный маршрут, чтобы убедиться в отсутствии подциклов:

function plot_tour(X, Y, x)
    plot = Plots.plot()
    for (i, j) in selected_edges(x, size(x, 1))
        Plots.plot!([X[i], X[j]], [Y[i], Y[j]]; legend = false)
    end
    return plot
end

plot_tour(X, Y, value.(iterative_model[:x]))

Метод отложенных ограничений

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

lazy_model = build_tsp_model(d, n)
function subtour_elimination_callback(cb_data)
    status = callback_node_status(cb_data, lazy_model)
    if status != MOI.CALLBACK_NODE_STATUS_INTEGER
        return  # Выполнить только при целочисленных решениях
    end
    cycle = subtour(callback_value.(cb_data, lazy_model[:x]))
    if !(1 < length(cycle) < n)
        return  # Добавлять ограничение только при наличии цикла
    end
    println("Found cycle of length $(length(cycle))")
    S = [(i, j) for (i, j) in Iterators.product(cycle, cycle) if i < j]
    con = @build_constraint(
        sum(lazy_model[:x][i, j] for (i, j) in S) <= length(cycle) - 1,
    )
    MOI.submit(lazy_model, MOI.LazyConstraint(cb_data), con)
    return
end
set_attribute(
    lazy_model,
    MOI.LazyConstraintCallback(),
    subtour_elimination_callback,
)
optimize!(lazy_model)
@assert is_solved_and_feasible(lazy_model)
objective_value(lazy_model)
525.7039004442727

Находится тот же оптимальный маршрут:

plot_tour(X, Y, value.(lazy_model[:x]))

Удивительно, но для этой конкретной модели с GLPK время решения хуже, чем при итеративном методе:

time_lazy = solve_time(lazy_model)
0.38207197189331055

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