Solves the traveling salesman problem using a callback for separating sub-tour elimination constraints.

The traveling salesman problem consists of finding a minimum-cost tour that visits each city exactly once and returns to the starting city.

Let \( V = \{1,\dots,n\} \) denote the set of cities, and let \( c_{ij} \) be the cost of traveling directly from city \( i \) to city \( j \).

An Exponential-Size Model

Decision Variables

We introduce binary variables

\[ x_{ij} = \begin{cases} 1 & \text{if the tour goes directly from city } i \text{ to city } j, \\ 0 & \text{otherwise}, \end{cases} \qquad \text{for all } i,j \in V,\; i \neq j. \]

Objective Function

The objective is to minimize the total travel cost of the tour:

\[ \min \sum_{i \in V} \sum_{j \in V \setminus \{i\}} c_{ij} x_{ij}. \]

Degree Constraints

To ensure that each city is entered and left exactly once, we impose the following constraints:

\[ \sum_{j \in V \setminus \{i\}} x_{ij} = 1 \quad \text{for all } i \in V, \]

\[ \sum_{i \in V \setminus \{j\}} x_{ij} = 1 \quad \text{for all } j \in V. \]

These constraints guarantee that every city has exactly one outgoing arc and one incoming arc. However, they still allow solutions consisting of multiple disconnected cycles (sub-tours).

Sub-Tour Elimination Constraints

To prevent sub-tours, we add the following family of constraints:

\[ \sum_{i \in S} \sum_{j \in S \setminus \{i\}} x_{ij} \le |S| - 1 \quad \text{for all } S \subset V,\; 2 \le |S| \le n-1. \]

Each constraint enforces that no proper subset \( S \) of cities can form a closed cycle. Together, these constraints ensure that the solution corresponds to a single Hamiltonian tour.

Since the number of subsets \( S \) grows exponentially with \( n \), these constraints are not added explicitly to the model. Instead, they are generated dynamically during the solution process.

Integrality Constraints

Finally, we require:

\[ x_{ij} \in \{0,1\} \quad \text{for all } i \neq j. \]

The Full Model

Finally, the full model is given by

\[ \begin{aligned} \min_x \quad & \sum_{i \in V} \sum_{j \in V \setminus \{i\}} c_{ij} x_{ij} \\ \text{s.t.}\quad & \sum_{j \in V \setminus \{i\}} x_{ij} = 1, \quad \text{for all } i \in V, \\ & \sum_{i \in V \setminus \{j\}} x_{ij} = 1, \quad \text{for all } j \in V, \\ & \sum_{i \in S} \sum_{j \in S \setminus \{i\}} x_{ij} \le |S| - 1, \quad \text{for all } S \subset V,\; 2 \le |S| \le n-1, \\ & x_{ij} \in \{0,1\}, \quad \text{for all } i \neq j. \end{aligned} \]

Implementation in idol

This example demonstrates how to solve the TSP using lazy constraints in Gurobi via idol. The model is initialized with the objective function and the degree constraints only. During the branch-and-bound process, a callback is used to inspect incumbent integer solutions. If a solution contains one or more sub-tours, the corresponding sub-tour elimination constraints are generated and added lazily.

//
// Created by henri on 20.12.25.
//
#include <iostream>
#include "idol/modeling.h"
#include "idol/mixed-integer/optimizers/branch-and-bound/BranchAndBound.h"
#include "idol/mixed-integer/optimizers/wrappers/Gurobi/Gurobi.h"
#include "idol/mixed-integer/optimizers/callbacks/Callback.h"
 
using namespace idol;
 
/*
 * This is the callback factory for the sub-tour elimination callback.
 * The role of this class is to create a callback object whenever it is needed.
 */
class SubTourEliminationCallback : public CallbackFactory {
    Vector<Var, 2> m_x; // We will need to have access to the x variables
public:
    SubTourEliminationCallback(const Vector<Var, 2>& t_x);
 
    Callback* operator()() override;
 
    [[nodiscard]] CallbackFactory* clone() const override;
 
    /*
     * This is the callback class.
     * The class overloads the operator() which is called whenever needed during the execution of the branch-and-bound method.
     */
    class Strategy;
};
 
class SubTourEliminationCallback::Strategy : public Callback {
    Vector<Var, 2> m_x; // original variables of the "master" problem
public:
    Strategy(const Vector<Var, 2>& t_x);
protected:
    void operator()(CallbackEvent t_event) override;
};
 
int main(int t_argc, const char** t_argv) {
 
    // City names (order matters)
    Vector<std::string> cities = {
        "Montpellier",
        "Grenoble",
        "Lille",
        "Pau",
        "Toulouse",
        "Paris",
        "Trier"
    };
 
    // Driving distance matrix in kilometers
    // matrix[i][j] = distance from cities[i] to cities[j]
    Vector<double, 2> distances = {
        {0.0,   296.0, 963.0, 780.0, 241.0, 750.0, 863.0}, // Montpellier
        {296.0,   0.0, 980.0, 860.0, 536.0, 570.0, 860.0}, // Grenoble
        {963.0, 980.0,   0.0,1000.0, 895.0, 225.0, 830.0}, // Lille
        {780.0, 860.0,1000.0,   0.0, 370.0, 795.0,1050.0}, // Pau
        {241.0, 536.0, 895.0, 370.0,   0.0, 680.0, 900.0}, // Toulouse
        {750.0, 570.0, 225.0, 795.0, 680.0,   0.0, 420.0}, // Paris
        {863.0, 860.0, 830.0,1050.0, 900.0, 420.0,   0.0}  // Trier
    };
 
    const auto n_cities = cities.size();
 
    /*********************************************/
    /* Create model without the lazy constraints */
    /*********************************************/
 
    Env env;
    Model model(env);
 
    const auto x = model.add_vars(Dim<2>(n_cities, n_cities), 0, 1, Binary, 0, "x");
 
    // Objective function
    model.set_obj_expr(idol_Sum(i, Range(n_cities), idol_Sum(j, Range(n_cities), distances[i][j] * x[i][j])));
 
    // Constraints
    for (unsigned int i = 0 ; i < n_cities ; ++i) {
        model.add_ctr(idol_Sum(j, Range(n_cities), x[i][j]) == 1);
        model.add_ctr(idol_Sum(j, Range(n_cities), x[j][i]) == 1);
    }
 
    // Remove x_ii since it is fixed to 0
    for (unsigned int i = 0 ; i < n_cities ; ++i) {
        model.remove(x[i][i]);
    }
 
    /*****************************************/
    /* Create optimizer and add the callback */
    /*****************************************/
 
    const auto gurobi = Gurobi()
                            .with_logs(true)
                            .with_lazy_cut(true) // this is necessary when using Gurobi with lazy constraints
                            .add_callback(SubTourEliminationCallback(x));
 
    model.use(gurobi);
 
    /*******************************/
    /* Optimize and print solution */
    /*******************************/
 
    model.optimize();
 
    for (unsigned int i = 0 ; i < n_cities ; ++i) {
        for (unsigned int j = 0 ; j < n_cities ; ++j) {
            if (i == j) { continue; }
            const bool is_one = model.get_var_primal(x[i][j]) > .5;
            if (is_one) {
                std::cout << cities[i] << " ---> " << cities[j] << std::endl;
            }
        }
    }
 
 
    return 0;
}
 
SubTourEliminationCallback::SubTourEliminationCallback(const Vector<Var, 2>& t_x) : m_x(t_x) {
 
}
 
CallbackFactory* SubTourEliminationCallback::clone() const {
    return new SubTourEliminationCallback(*this);
}
 
Callback* SubTourEliminationCallback::operator()() {
    return new Strategy(m_x);
}
 
SubTourEliminationCallback::Strategy::Strategy(const Vector<Var, 2>& t_x) : m_x(t_x) {
 
}
 
void SubTourEliminationCallback::Strategy::operator()(CallbackEvent t_event) {
 
    // If the event that triggered the callback is not a new incumbent solution being found, we directly exit the callback.
    if (t_event != IncumbentSolution) {
        return;
    }
 
    const unsigned int n_cities = m_x.size();
    std::vector<int> visited(n_cities, 0);
 
    // Find all cycles in the incumbent
    const auto incumbent = primal_solution();
    for (int start = 0; start < n_cities; ++start) {
 
        if (visited[start]) continue;
 
        std::vector<int> cycle;
        int curr = start;
 
        while (!visited[curr]) {
            visited[curr] = 1;
            cycle.push_back(curr);
 
            // Move to the next city in the incumbent solution
            int next = -1;
            for (int j = 0; j < n_cities; ++j) {
                if (curr == j) { continue; }
                if (incumbent.get(m_x[curr][j]) > 0.5) {
                    next = j;
                    break;
                }
            }
 
            if (next == -1) break; // isolated node (should not happen)
            curr = next;
        }
 
        // Add a subtour elimination cut if this is a proper subtour
        if (cycle.size() < n_cities) {
            LinExpr cut_lhs;
            for (int i : cycle) {
                for (int j : cycle) {
                    if (i != j) {
                        cut_lhs += m_x[i][j];
                    }
                }
            }
 
            std::cout << "Adding subtour cut for " << cycle.size() << " cities" << std::endl;
            add_lazy_cut(cut_lhs <= cycle.size() - 1);
 
        }
    }
 
}