[PYTHON] About the Ordered Traveling Salesman Problem

Have you ever heard of the traveling salesman problem? The problem is to find the circuit that minimizes the path length among the circuits that can trace the graph with a single stroke.

This time, I would like to deal with the special ** Sequential Ordering Problem ** among the traveling salesman problems. (I can't find a Japanese translation, so I named it appropriately, so if anyone knows it, please let me know.)

What is different from the original traveling salesman problem is that there is an order constraint that "this point must go after this point". This makes it possible to consider situations such as "giving the baggage received from Mr. A to Mr. B".

This time, we will deal with the formulation and numerical experiment results using the integer programming of the Sequential Ordering Problem.

First, prepare symbols and constants.

V :Point set, |V| := n \\
ss \in V :start point\\ 
E :Branch set\\
c_{ij} ((v_i,v_j)\in E) :branch(v_i,v_j)Cost to go through\\
\mathcal{S} = \{[v_i,v_j], \cdots\} :A set with a list of two points that must be in order

Next, define the variables.

x_{ij} = 
\begin{cases}
    1 &branch(v_i,v_j)When passing through\\
    0 & otherwise
  \end{cases}\\
1\le u_{i} \le n-1 :Point v_Order through i

Now let's move on to the formulation. The meaning of each detailed constraint will be added later.

\begin{eqnarray}
{minimize}& \sum_{(v_i,v_j)\in E}c_{ij}x_{ij} \\
subject \ to & \ \sum_{i} x_{ij} = 1 (\forall j) \\
&\sum_{j} x_{ij} = 1 (\forall i)\\
&u_i - u_j + (n-1)x_{ij} \le n-2 \ (\forall j \in S \setminus \{ss\}) \\
&u_{ss} = 0 \\
&u_{i} + 1 \le u_j \ (\forall (v_i,v_j) \in \mathcal{S})
\end{eqnarray}

This is an implementation example. I made a random graph and experimented.

sop.py

import pulp
import random
import networkx as nx
import math,time
import matplotlib.pyplot as plt

def make_random_graph(n):
    G = nx.DiGraph()
    for i in range(n):
        G.add_node(i,x=random.randint(0,10),y=random.randint(0,10))
    for i in range(n):
        for j in range(n):
            if i != j:
                dist = math.sqrt((G.nodes()[i]['x']-G.nodes()[j]['x'])**2 + (G.nodes()[i]['y']-G.nodes()[j]['y'])**2)
                G.add_edge(i,j,dist=dist)
    return G

def get_random_sequential_order(num_node,m):
    box = set()
    while len(box) < m:
        i = random.randint(0,num_node-2)
        j = random.randint(i+1,num_node-1)
        if (i,j) not in box:
            box.add((i,j))
    return box

def solve_SOP(G,precedense,num_node,ss):
    problem = pulp.LpProblem(name='SOP',sense=pulp.LpMinimize)
    x = {(i,j):pulp.LpVariable(cat="Binary",name=f"x_{i}_{j}") for (i,j) in G.edges()}
    u = {i:pulp.LpVariable(cat="Integer",name=f"u_{i}",lowBound=1,upBound=num_node-1) for i in G.nodes()}
    cost = {(i,j):G.adj[i][j]['dist'] for (i,j) in G.edges()}

    problem += pulp.lpSum([x[(i,j)]*cost[(i,j)] for (i,j) in G.edges()])


    for i in G.nodes():
        problem.addConstraint(pulp.lpSum([x[(i,j)] for j in range(num_node) if j != i]) == 1, f'outflow_{i}')
        problem.addConstraint(pulp.lpSum([x[(j,i)] for j in range(num_node) if j != i]) == 1, f'inflow_{i}')

    for i,j in G.edges():
        if i != ss and j != ss:
            problem.addConstraint(u[i]-u[j]+(num_node-1)*x[i,j] <= num_node-2, f'up_{i}_{j}')

    for i,j in precedense:
        problem.addConstraint(u[i]+1 <= u[j], f'sequential_{i}_{j}')

    u[ss] = 0

    print('start solving')
    start = time.time()
    status = problem.solve(pulp.CPLEX())
    # status = problem.solve()
    print(pulp.LpStatus[status])

    print(time.time()-start)

    if pulp.LpStatus[status] != 'Optimal':
        print('Infeasible!')
        exit()

    return x,u

def plot(G,x,u,precedense,ss):
    pos = {i: (G.nodes()[i]['x'], G.nodes()[i]['y']) for i in G.nodes()}
    nx.draw_networkx_nodes(G, pos, node_size=100, alpha=1, node_color='skyblue')
    edgelist = [e for e in G.edges() if x[e].value() > 0]
    nx.draw_networkx_edges(G, pos, edgelist=edgelist,width=3)
    precedense = [e for e in precedense]
    nx.draw_networkx_edges(G, pos, edgelist=precedense,edge_color='red')
    for i in G.nodes():
        if i != ss:
            plt.text(G.nodes()[i]['x'],G.nodes()[i]['y'],int(u[i].value()))
        else:
            plt.text(G.nodes()[i]['x'],G.nodes()[i]['y'],u[i])


    plt.show()

def main():
    num_node = 10
    num_precedence = 5
    ss = 0
    precedense = get_random_sequential_order(num_node,num_precedence)
    print(precedense)
    G = make_random_graph(num_node)
    x,u = solve_SOP(G,precedense,num_node,ss)
    plot(G,x,u,precedense,ss)

if __name__ == '__main__':
    main()

It is the execution result with 10 points.

The black arrow indicates the order of the circuit, and the red arrow indicates the preference order. As the number of points increases, the calculation time increases exponentially, so it seems to be difficult for the branch-and-bound method that finds the exact solution of the integer programming method. Therefore, an approximate solution has been proposed as well as the original traveling salesman problem.