tsp_ilp.py

import math
from builtins import min

import scipy.spatial as sp
import pandas as pd
import numpy as np
import logging as log
import matplotlib.pyplot as plt
import json
import sys
# also required packages: openpyxl
from ortools.linear_solver import pywraplp


def solve_OrTools(dima):
    '''
    Generate MIP (Mixed-Integer Programming) model using Google OR-Tools and solve it.

    :param dima: the distance matrix
    :return: solution X, model, status
    '''

    if dima.ndim != 2 or dima.shape[0] != dima.shape[1]:
        raise ValueError("Invalid dima dimensions detected. Square matrix expected.")

    # Determine the number of nodes
    num_nodes = dima.shape[0]
    all_nodes = range(0, num_nodes)
    all_but_first_nodes = range(1, num_nodes)

    # Create the model.
    solver_name = 'CP-SAT'
    log.info('Instantiating solver ' + solver_name)
    model = pywraplp.Solver.CreateSolver(solver_name)
    model.EnableOutput()

    log.info('Defining MIP model... ')
    # Generating decision variables X_ij
    log.info('Creating ' + str(num_nodes * num_nodes) + ' boolean x_ij variables... ')
    x = {}
    for i in all_nodes:
        for j in all_nodes:
            x[(i, j)] = model.BoolVar('x_i%ij%i' % (i, j))

    log.info('Creating ' + str(num_nodes) + ' boolean u_i variables... ')
    u = {}
    for i in all_nodes:
        u[i] = model.IntVar(0, num_nodes, 'u_i%i' % i)

    # Constraint 1: Leave every point exactly once
    log.info('Creating ' + str(num_nodes) + ' Constraint 1... ')
    for i in all_nodes:
        model.Add(sum(x[(i, j)] for j in all_nodes) == 1)

    # Constraint 2: Reach every point from exactly one other point
    log.info('Creating ' + str(num_nodes) + ' Constraint 2... ')
    for j in all_nodes:
        model.Add(sum(x[(i, j)] for i in all_nodes) == 1)

    # Constraint 3.1: Subtour elimination constraints (Miller-Tucker-Zemlin) part 1
    log.info('Creating 1 Constraint 3.1... ')
    model.Add(u[0] == 1)

    # Constraint 3.2: Subtour elimination constraints (Miller-Tucker-Zemlin) part 2
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        model.Add(2 <= u[i])
        model.Add(u[i] <= num_nodes)

    # Constraint 3.3: Subtour elimination constraints (Miller-Tucker-Zemlin) part 3
    log.info('Creating ' + str(len(all_but_first_nodes)) + ' Constraint 3.2... ')
    for i in all_but_first_nodes:
        for j in all_but_first_nodes:
            model.Add(u[i] - u[j] + 1 <= (num_nodes - 1) * (1 - x[(i, j)]))

    # Minimize the total distance
    model.Minimize(sum(x[(i, j)] * dima[(i, j)] for i in all_nodes for j in all_nodes))

    log.info('Solving MIP model... ')
    status = model.Solve()

    return x, model, status


def print_solution(x):
    num_nodes = len(x)
    all_nodes = range(0, int(math.sqrt(num_nodes)))
    for i in all_nodes:
        for j in all_nodes:
            if int(x[(i, j)].solution_value()) == 1:
                log.info(f'x({i}/{j})=' + str(int(x[i, j].solution_value())))


def extract_tsp_sequence(x, num_nodes):
    """Extracts the order of nodes from the decision variables x[i,j].

    Args:
        x: A dictionary with the 2D decision variables x[i,j].
        num_nodes: The number of nodes in the TSP.

    Returns:
        A list containing the order of visited nodes.
    """
    # Reverse mapping for the next node: from node to node
    next_node = {}

    # Fill the next_node map based on the decision variables
    for i in range(num_nodes):
        for j in range(num_nodes):
            if x[i, j].solution_value() == 1:
                next_node[i] = j

    # Find the order of visited nodes
    sequence = []
    current_node = 0  # We start with a starting node, e.g., 0
    for _ in range(num_nodes):
        sequence.append(current_node)
        current_node = next_node[current_node]

    return sequence


def plot_tsp_solution(tsp, solution, title):
    # Extract coordinates
    x_coords = [tsp.lat[node] for node in solution]
    y_coords = [tsp.lng[node] for node in solution]

    # Add the start point to the end to complete the tour
    x_coords.append(x_coords[0])
    y_coords.append(y_coords[0])

    # Plot the nodes
    plt.scatter(x_coords, y_coords, c='gray')

    # Plot the path
    plt.plot(x_coords, y_coords, c='green')

    # Annotate the nodes
    if len(tsp) < 200:
        for idx, node in enumerate(solution):
            plt.annotate(str(node), (x_coords[idx], y_coords[idx]))

    plt.xlabel('X Coordinate')
    plt.ylabel('Y Coordinate')
    plt.title('TSP Visualization (' + title + ')')
    plt.show()


def main():
    # Configure logger for info level
    log.basicConfig(
        format='%(asctime)s %(levelname)-8s %(message)s',
        level=log.INFO,
        datefmt='%Y-%m-%d %H:%M:%S',
        stream=sys.stdout)

    # Load TSP instance
    tsp_problem = 'berlin52.tsp'

    log.info("Reading TSP problem Instance " + tsp_problem)
    tsp = pd.read_csv('./data/' + tsp_problem, sep=' ', skiprows=6, dtype=float,
                      names=['nodeId', 'lat', 'lng'], skipfooter=1, engine='python')
    tsp = tsp.sort_values(by='nodeId', inplace=False)

    A = tsp[['lat', 'lng']].to_numpy()
    dima = sp.distance_matrix(A, A)

    # Now solve the problem
    x, model, status = solve_OrTools(dima)

    # Check problem response
    if status == pywraplp.Solver.OPTIMAL:
        log.info('Solution:')
        log.info('Optimal solution found.')
        log.info('Objective value =' + str(model.Objective().Value()))
        tour = extract_tsp_sequence(x, len(A))
        print(f"Optimal Tour: {tour}")
        plot_tsp_solution(tsp, tour, "berlin52")
    elif status == pywraplp.Solver.INFEASIBLE:
        log.info('The problem is infeasible.')
    else:
        log.info('The problem could not be solved. Return state was: ' + status)


# Press the green button in the gutter to run the script.
if __name__ == '__main__':
    main()