時間枠付き巡回セールスマン問題に対する定式化とアルゴリズム

準備

from gurobipy import Model, quicksum, GRB
import math
import networkx as nx

時間枠付き巡回セールスマン問題

ここでは,巡回セールスマン問題に時間枠を追加した 時間枠付き巡回セールスマン問題(traveling salesman problem with time windows)を考える.

この問題は,特定の点 $0$ を時刻 $0$ に出発すると仮定し, 点間の移動距離 $c_{ij}$ を移動時間とみなし, さらに点 $i$ に対する出発時刻が最早時刻 $e_i$ と最遅時刻 $\ell_i$ の間でなければならないという制約を課した問題である. ただし,時刻 $e_i$ より早く点 $i$ に到着した場合には,点 $i$ 上で時刻 $e_i$ まで待つことができるものとする.

ポテンシャル定式化

巡回セールスマン問題に対するポテンシャル制約の拡張を考える.

点 $i$ を出発する時刻を表す変数 $t_i$ を導入する. $t_i$ は以下の制約を満たす必要がある. $$ e_i \leq t_i \leq \ell_i \ \ \ \forall i=1,2,\ldots,n $$ ただし, $e_1=0, \ell_1=\infty$ と仮定する.

点 $i$ の次に点 $j$ を訪問する $(x_{ij}=1)$ ときには, 点 $j$ を出発する時刻 $t_j$ は,点 $i$ を出発する時刻に移動時間 $c_{ij}$ を加えた値以上であることから, 以下の式を得る. $$  t_i + c_{ij} - M (1-x_{ij}) \leq t_j \ \ \ \forall i,j, j \neq 1, i \neq j $$ ここで,$M$ は大きな数を表す定数である. なお,移動時間 $c_{ij}$ は正の数と仮定する.$c_{ij}$ が $0$ だと $t_i=t_j$ になる可能性があり, 部分巡回路ができてしまう.これを避けるためには,巡回セールスマン問題と同様の制約を付加する必要があるが, $c_{ij}>0$ の仮定の下では,上の制約によって部分巡回路を除去することができる.

このような大きな数Big Mを含んだ定式化はあまり実用的ではないので,時間枠を用いて強化したものを示す.

$$ \begin{array}{lll} minimize & \sum_{i \neq j} c_{ij} x_{ij} & \\ s.t. & \sum_{j: j \neq i} x_{ij} = 1 & \forall i=1,2,\ldots,n \\ & \sum_{j: j \neq i} x_{ji} = 1 & \forall i=1,2,\ldots,n \\ & t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) \leq t_j & \forall i,j,j \neq 1, i \neq j \\ & x_{ij} \in \{0,1\} & \forall i,j, i \neq j \\ & e_i \leq t_{i} \leq \ell_i & \forall i=1,2,\ldots,n \end{array} $$

巡回セールスマン問題のときと同様に,ポテンシャル制約と上下限制約は, 持ち上げ操作によってさらに以下のように強化できる.

$$ t_i + c_{ij} - [\ell_i +c_{ij}-e_j]^+ (1-x_{ij}) + [\ell_i-e_j+\min \{-c_{ji}, e_j-e_i\} ]^+ x_{ji} \leq t_j \ \ \ \forall i,j, j \neq1, i \neq j $$$$ e_i + \sum_{j \neq i} [e_j+ c_{ji} -e_i]^{+} x_{ji} \leq t_i \ \ \ \forall i=2,\cdots,n $$$$ t_i \leq \ell_i -\sum_{j \neq 1,i} [\ell_i -\ell_j+c_{ij}]^+ x_{ij} \ \ \ \forall i=2,\cdots,n $$

以下に,強化していない標準定式化 mtztw と強化した定式化 mtz2tw のコードを示す.

def mtztw(n, c, e, l):
    """mtzts: model for the traveling salesman problem with time windows
    (based on Miller-Tucker-Zemlin's one-index potential formulation)
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
        - e[i]: earliest date for visiting node i
        - l[i]: latest date for visiting node i
    Returns a model, ready to be solved.
    """
    model = Model("tsptw - mtz")
    x, u, v = {}, {}, {}
    for i in range(1, n + 1):
        u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
        v[i] = model.addVar(lb=0, ub=n - 1, vtype="C", name="v(%s)" % i)
        for j in range(1, n + 1):
            if i != j:
                x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
    model.update()

    for i in range(1, n + 1):
        model.addConstr(
            quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
        )
        model.addConstr(
            quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
        )

    for i in range(1, n + 1):
        for j in range(2, n + 1):
            if i != j:
                M = max(l[i] + c[i, j] - e[j], 0)
                model.addConstr(
                    u[i] - u[j] + M * x[i, j] <= M - c[i, j], "MTZ(%s,%s)" % (i, j)
                )
                model.addConstr(
                    v[i] - v[j] + (n - 1) * x[i, j] <= n - 2, "MTZorig(%s,%s)" % (i, j)
                )

    model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)

    model.update()
    model.__data = x, u
    return model


def mtz2tw(n, c, e, l):
    """mtz: model for the traveling salesman problem with time windows
    (based on Miller-Tucker-Zemlin's one-index potential formulation, stronger constraints)
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
        - e[i]: earliest date for visiting node i
        - l[i]: latest date for visiting node i
    Returns a model, ready to be solved.
    """
    model = Model("tsptw - mtz-strong")
    x, u = {}, {}
    for i in range(1, n + 1):
        u[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="u(%s)" % i)
        for j in range(1, n + 1):
            if i != j:
                x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
    model.update()

    for i in range(1, n + 1):
        model.addConstr(
            quicksum(x[i, j] for j in range(1, n + 1) if j != i) == 1, "Out(%s)" % i
        )
        model.addConstr(
            quicksum(x[j, i] for j in range(1, n + 1) if j != i) == 1, "In(%s)" % i
        )

        for j in range(2, n + 1):
            if i != j:
                M1 = max(l[i] + c[i, j] - e[j], 0)
                M2 = max(l[i] + min(-c[j, i], e[j] - e[i]) - e[j], 0)
                model.addConstr(
                    u[i] + c[i, j] - M1 * (1 - x[i, j]) + M2 * x[j, i] <= u[j],
                    "LiftedMTZ(%s,%s)" % (i, j),
                )

    for i in range(2, n + 1):
        model.addConstr(
            e[i]
            + quicksum(
                max(e[j] + c[j, i] - e[i], 0) * x[j, i]
                for j in range(1, n + 1)
                if i != j
            )
            <= u[i],
            "LiftedLB(%s)" % i,
        )

        model.addConstr(
            u[i]
            <= l[i]
            - quicksum(
                max(l[i] - l[j] + c[i, j], 0) * x[i, j]
                for j in range(2, n + 1)
                if i != j
            ),
            "LiftedUB(%s)" % i,
        )

    model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)

    model.update()
    model.__data = x, u
    return model

ベンチマーク問題例の読み込み

時間枠付き巡回セールスマン問題のベンチマーク問題例は,以下のサイトから入手できる.

https://lopez-ibanez.eu/tsptw-instances

読み込んでデータを準備する.

folder = "../data/tsptw/"
# fn = "n60w20.001.txt"
fn ="n200w20.002.txt"
# fn ="n100w40.002.txt"
# fn ="n100w60.002.txt"
f = open(folder + fn)
data = f.readlines()
f.close()
x_, y_, early, late = [], [], [], []
for row in data:
    try:
        L = list(map(float, row.split()))
        if L[0] >= 900:
            continue
        x_.append(L[1])
        y_.append(L[2])
        early.append(L[4])
        late.append(L[5])
    except:
        pass
n = len(x_)
e, l = {}, {}
for i in range(n):
    e[i + 1] = int(early[i])
    l[i + 1] = int(late[i])
c = {}
for i in range(n):
    for j in range(n):
        c[i + 1, j + 1] = int(math.sqrt((x_[i] - x_[j]) ** 2 + (y_[i] - y_[j]) ** 2))
for i in range(1, n + 1):
    for j in range(1, n + 1):
        for k in range(1, n + 1):
            c[i, j] = min(c[i, j], c[i, k] + c[k, j])

求解と描画

まず,強化していない標準定式化を用いて求解する.

model = mtztw(n, c, e, l)
model.optimize()
x, u = model.__data

sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]

print("Opt.value =", model.ObjVal, sol)
Academic license - for non-commercial use only - expires 2024-11-13
Using license file /Users/mikiokubo/gurobi.lic
Gurobi Optimizer version 9.1.1 build v9.1.1rc0 (mac64)
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 7322 rows, 3782 columns and 27369 nonzeros
Model fingerprint: 0x6be2b77d
Variable types: 122 continuous, 3660 integer (3660 binary)
Coefficient statistics:
  Matrix range     [1e+00, 6e+02]
  Objective range  [1e+00, 6e+01]
  Bounds range     [1e+00, 6e+02]
  RHS range        [1e+00, 6e+02]
Presolve removed 4550 rows and 2423 columns
Presolve time: 0.11s
Presolved: 2772 rows, 1359 columns, 18386 nonzeros
Variable types: 112 continuous, 1247 integer (1247 binary)

Root relaxation: objective 5.312981e+02, 313 iterations, 0.01 seconds

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  531.29814    0   40          -  531.29814      -     -    0s
     0     0  539.83333    0   35          -  539.83333      -     -    0s
     0     0  539.83333    0   35          -  539.83333      -     -    0s
     0     0  541.00000    0   27          -  541.00000      -     -    0s
     0     0  541.00000    0   30          -  541.00000      -     -    0s
H    0     0                     555.0000000  541.00000  2.52%     -    0s
H    0     0                     554.0000000  541.00000  2.35%     -    0s
     0     0  541.00000    0    9  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   38  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   25  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   26  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   24  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   24  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0    6  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0   10  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0    6  554.00000  541.00000  2.35%     -    0s
     0     0  541.00000    0    6  554.00000  541.00000  2.35%     -    0s
     0     0  541.25000    0    6  554.00000  541.25000  2.30%     -    0s
     0     0  541.25000    0    6  554.00000  541.25000  2.30%     -    0s
     0     0  541.27778    0   19  554.00000  541.27778  2.30%     -    0s
     0     0  541.27778    0   17  554.00000  541.27778  2.30%     -    0s
     0     0  541.27778    0   21  554.00000  541.27778  2.30%     -    0s
     0     0  541.27778    0   17  554.00000  541.27778  2.30%     -    0s
     0     0  542.48649    0    9  554.00000  542.48649  2.08%     -    0s
     0     0  542.67925    0    9  554.00000  542.67925  2.04%     -    0s
     0     0  543.07143    0    8  554.00000  543.07143  1.97%     -    0s
H    0     0                     552.0000000  543.07143  1.62%     -    0s
     0     0  543.13514    0   10  552.00000  543.13514  1.61%     -    0s
     0     0  543.63975    0    9  552.00000  543.63975  1.51%     -    0s
     0     0  543.75000    0   10  552.00000  543.75000  1.49%     -    0s
     0     0  544.00000    0    8  552.00000  544.00000  1.45%     -    1s
     0     0  544.11951    0   10  552.00000  544.11951  1.43%     -    1s
H    0     0                     551.0000000  544.39024  1.20%     -    1s
     0     0  544.39024    0    6  551.00000  544.39024  1.20%     -    1s
     0     0  544.42857    0   10  551.00000  544.42857  1.19%     -    1s
     0     0  544.68750    0    9  551.00000  544.68750  1.15%     -    1s
     0     0  544.68750    0    9  551.00000  544.68750  1.15%     -    1s
     0     2  544.68750    0    9  551.00000  544.68750  1.15%     -    1s

Cutting planes:
  Learned: 5
  Gomory: 3
  Cover: 1
  Implied bound: 16
  Clique: 1
  MIR: 19
  StrongCG: 2
  Inf proof: 3
  RLT: 1
  Relax-and-lift: 13

Explored 329 nodes (4493 simplex iterations) in 1.20 seconds
Thread count was 16 (of 16 available processors)

Solution count 4: 551 552 554 555 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.510000000000e+02, best bound 5.510000000000e+02, gap 0.0000%
Opt.value = 551.0 [1, 7, 13, 39, 49, 10, 16, 52, 35, 57, 18, 20, 43, 59, 23, 60, 33, 8, 26, 14, 5, 53, 6, 51, 45, 31, 48, 19, 47, 32, 17, 24, 22, 3, 29, 42, 11, 30, 12, 28, 40, 44, 46, 36, 27, 9, 58, 4, 34, 55, 2, 25, 50, 38, 21, 37, 56, 15, 41, 61, 54]

次に,時間枠を用いて制約を強化した定式化を用いて求解する.

model = mtz2tw(n, c, e, l)
model.optimize()
x, u = model.__data

sol = [i for (v, i) in sorted([(u[i].X, i) for i in u])]
print("Opt.value =", model.ObjVal, sol)
Gurobi Optimizer version 9.5.2 build v9.5.2rc0 (mac64[x86])
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 3842 rows, 3721 columns and 22227 nonzeros
Model fingerprint: 0x10540bf3
Variable types: 61 continuous, 3660 integer (3660 binary)
Coefficient statistics:
  Matrix range     [1e+00, 6e+02]
  Objective range  [1e+00, 6e+01]
  Bounds range     [1e+00, 6e+02]
  RHS range        [1e+00, 6e+02]
Presolve removed 2118 rows and 2215 columns
Presolve time: 0.05s
Presolved: 1724 rows, 1506 columns, 7938 nonzeros
Variable types: 57 continuous, 1449 integer (1449 binary)

Root relaxation: objective 5.361061e+02, 336 iterations, 0.01 seconds (0.01 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

     0     0  536.10606    0   67          -  536.10606      -     -    0s
     0     0  541.00000    0   13          -  541.00000      -     -    0s
H    0     0                     553.0000000  541.00000  2.17%     -    0s
     0     0  541.00000    0   10  553.00000  541.00000  2.17%     -    0s
H    0     0                     552.0000000  541.00000  1.99%     -    0s
H    0     0                     551.0000000  541.00000  1.81%     -    0s
     0     0  541.00000    0    6  551.00000  541.00000  1.81%     -    0s
     0     0  541.00000    0   42  551.00000  541.00000  1.81%     -    0s
     0     0  541.00000    0   16  551.00000  541.00000  1.81%     -    0s
     0     0  541.05935    0   11  551.00000  541.05935  1.80%     -    0s
     0     0  541.05935    0   11  551.00000  541.05935  1.80%     -    0s
     0     0  541.37202    0   11  551.00000  541.37202  1.75%     -    0s
     0     0  541.42378    0   11  551.00000  541.42378  1.74%     -    0s
     0     0  541.76983    0   16  551.00000  541.76983  1.68%     -    0s
     0     0  541.78378    0   18  551.00000  541.78378  1.67%     -    0s
     0     0  542.14093    0   18  551.00000  542.14093  1.61%     -    0s
     0     0  542.21128    0   17  551.00000  542.21128  1.60%     -    0s
     0     0  542.22222    0   26  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   11  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   13  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   11  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   11  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   11  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   11  551.00000  542.22222  1.59%     -    0s
     0     0  542.22222    0   12  551.00000  542.22222  1.59%     -    0s
     0     0  542.50824    0   11  551.00000  542.50824  1.54%     -    0s
     0     0  542.50824    0   10  551.00000  542.50824  1.54%     -    0s
     0     2  542.50824    0   10  551.00000  542.50824  1.54%     -    0s

Cutting planes:
  Learned: 5
  Gomory: 5
  Cover: 2
  MIR: 4
  StrongCG: 2
  GUB cover: 2
  Zero half: 5
  RLT: 3
  Relax-and-lift: 3

Explored 348 nodes (3140 simplex iterations) in 0.48 seconds (0.26 work units)
Thread count was 16 (of 16 available processors)

Solution count 3: 551 552 553 

Optimal solution found (tolerance 1.00e-04)
Best objective 5.510000000000e+02, best bound 5.510000000000e+02, gap 0.0000%
Opt.value = 551.0 [1, 7, 13, 39, 49, 10, 16, 52, 35, 57, 18, 20, 43, 59, 23, 33, 60, 8, 26, 14, 5, 53, 6, 45, 51, 31, 48, 19, 47, 32, 17, 24, 22, 3, 29, 42, 11, 30, 12, 28, 40, 44, 46, 36, 27, 9, 58, 4, 34, 55, 2, 25, 50, 38, 21, 37, 56, 15, 41, 61, 54]
G = nx.Graph()
for idx, i in enumerate(sol[:-1]):
    G.add_edge(i, sol[idx + 1])
G.add_edge(sol[-1],sol[1])
pos = {i: (x_[i-1], y_[i-1]) for i in range(1,n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue");