巡回セールスマン問題について耳にしたことはございますか？
グラフを一筆書きで辿れる巡回路の中で、経路長が最小になる巡回路を求める問題です。

今回はその巡回セールスマン問題の中でも特殊な**順序付き巡回セールスマン問題(Sequential Ordering Problem)**について扱いたいと思います。(日本語訳が見つからなくて適当に名付けているのでご存知の方がいれば教えて頂けると助かります。)

元の巡回セールスマン問題と何が違うのかというと、「この点はこの点よりも後に行かなければならない」という順序に関する制約がついていることです。
これによって例えば、「Aさんから預かった荷物をBさんに渡す」といった状況を考慮することができます。

今回はそのSequential Ordering Problemの整数計画を用いた定式化及び数値実験結果について扱います。

まずは記号及び定数の準備をします。

V : 点集合, |V| := n \\
ss \in V : 始点\\ 
E : 枝集合 \\
c_{ij} ((v_i,v_j)\in E) : 枝(v_i,v_j)を通るコスト\\
\mathcal{S} = \{[v_i,v_j], \cdots\} : 順序を守らなければならない２点のリストを持つ集合

次に変数を定義します。

x_{ij} = 
\begin{cases}
    1 & 枝(v_i,v_j)を通るとき \\
    0 & otherwise
  \end{cases}\\
1\le u_{i} \le n-1 : 点v_iを通る順番

では定式化に移ります。細かいそれぞれの制約の意味は後日追記します。

\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 V \setminus \{ss\}) \\
&u_{ss} = 0 \\
&u_{i} + 1 \le u_j \ (\forall (v_i,v_j) \in \mathcal{S})
\end{eqnarray}

実装例です。
ランダムなグラフを作って実験してみました。

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()

10点での実行結果です。

黒矢印が求まった巡回路の順番で、赤矢印が選好順序を表しています。
点数を増やすと計算時間が指数的に大きくなるので、整数計画法の厳密解を求める分枝限定法では厳しそうです。それ故、オリジナルの巡回セールスマン問題と同様近似解法が提案されています。