効率的な配送を実現せよ：Pythonで学ぶ配送ルート最適化

物流コストの削減と顧客満足度の向上は、多くの企業にとって重要な課題です。

今回は、Pythonを使った数理最適化モデルを構築し、最適な配送ルートを決定する方法をお話しします。

巡回セールスマン問題（TSP）から始め、より複雑な車両経路問題（VRP）まで、段階的にお話ししていきます。

この手法を身につければ、データに基づいた効率的な配送計画の立案が可能になるでしょう。

はじめに

配送ルートの最適化は、物流業界だけでなく、あらゆる配送業務を行う企業にとって重要な課題です。

適切なルート選択により、燃料コストの削減、配送時間の短縮、車両の効率的な利用が可能になります。しかし、配送先の増加や時間枠制約、車両の容量制限など、考慮すべき要素は多岐にわたります。

ここで数理最適化の出番です。数理最適化を用いることで、複雑な制約条件を考慮しつつ、最適な配送ルートを導き出すことができます。

基本的な巡回セールスマン問題（TSP）から始め、より現実的な車両経路問題（VRP）へと発展させていきます。

それぞれの問題について、以下の要素を説明します。

問題設定
数理モデル
Pythonによる実装

最適化ライブラリの概要：PuLPとOR-Tools

今回は以下の2つのライブラリを使います。

PuLP：比較的シンプルな問題
OR-Tools：より複雑で大規模な問題

両ライブラリとも強力な最適化ツールですが、問題の特性や要件に応じて適切なものを選択することが重要です。

　PuLP

PuLPは、Pythonで線形計画問題（LP）と整数計画問題（IP）を解くためのオープンソースライブラリです。

pipでインストールするときのコードは以下です。

pip install pulp

pip install pulp

特徴

使いやすいPythonインターフェース
多様なソルバーをサポート（CBC, GLPK, CPLEX, Gurobi等）
問題の定式化が直感的
中小規模の問題に適している

使用例

巡回セールスマン問題（TSP）の解法
生産計画最適化
資源配分問題

　OR-Tools

OR-Toolsは、最適化問題を解くためのオープンソースソフトウェアスイートです。

pipでインストールするときのコードは以下です。

pip install ortools

pip install ortools

特徴

線形計画、整数計画、制約充足問題、ネットワークフロー等、多様な最適化問題に対応
高性能なソルバーを内蔵
大規模な問題にも対応可能
複雑な制約を扱うのに適している

使用例

車両経路問題（VRP）の解法
スケジューリング問題
配送最適化
施設配置問題

　2つのライブラリの選択基準

問題の種類

単純な線形計画問題 → PuLP
複雑な制約を持つ問題や組合せ最適化問題 → OR-Tools

問題の規模

小〜中規模の問題 → PuLP
大規模な問題 → OR-Tools

使いやすさ

初心者や簡単な問題の場合 → PuLP
より高度な最適化や制御が必要な場合 → OR-Tools

実行速度

一般的に、OR-Toolsの方が高速

多くの場合、PuLPで始めて、より複雑な問題や大規模な問題に直面した際にOR-Toolsに移行するのが良いです。

巡回セールスマン問題（TSP）

　問題設定

ある配送会社の事例です。

配送車ごとに、その日の配送先である顧客が決められています。

最短ルートで顧客を回ろうと考えています。

与えられた情報

配送拠点（デポ）：1箇所
顧客数：n
各地点間の距離が既知

目的

全ての顧客を1回ずつ訪問し、デポに戻ってくる最短ルートを見つける

決定変数

$x_{ij}$ : 地点 i から地点 j へ移動する場合は 1、そうでない場合は 0

制約条件

各顧客を1回だけ訪問すること
デポから出発し、デポに戻ってくること
部分巡回（サブツアー）の除去

　数理モデル

　　モデル全体

巡回セールスマン問題（TSP: Traveling Salesman Problem）を定式化すると、以下のようになります。

\begin{array}{ll}\displaystyle \text{Minimize} \quad & \displaystyle \sum_{i=0}^{n} \sum_{j=0, j \neq i}^{n} d_{ij} x_{ij} \\\text{subject to} \quad & \displaystyle \sum_{j=0, j \neq i}^{n} x_{ij} = 1 \quad \forall i \in \{0, 1, \ldots, n\} \\& \displaystyle \sum_{j=0, j \neq i}^{n} x_{ji} = 1 \quad \forall i \in \{0, 1, \ldots, n\} \\[0.5cm] & \displaystyle u_i - u_j + n x_{ij} \leq n - 1 \quad \forall i, j \in \{1, 2, \ldots, n\}, i \neq j \\[0.3cm] & \displaystyle x_{ij} \in \{0, 1\} \quad \forall i, j \in \{0, 1, \ldots, n\}, i \neq j \\[0.3cm] & \displaystyle u_i \in \{0, 1, \ldots, n\} \quad \forall i \in \{1, 2, \ldots, n\}\end{array}

以下、変数です。

$x_{ij} \in \{0, 1\}$ : 地点 $i$ から地点 $j$ への移動が選択された場合に1、そうでない場合に0。
$u_i \in \{0, 1, \ldots, n\}$ : 地点 $i$ の順序を表す補助変数。

　　目的関数

総移動距離を最小化することを目的としています。

\displaystyle \text{Minimize} \quad \sum_{i=0}^{n} \sum_{j=0, j \neq i}^{n} d_{ij} x_{ij}

ここで、 $d_{ij}$ は地点 $i$ から地点 $j$ への距離を表し、 $x_{ij}$ は地点 $i$ から地点 $j$ への移動が選択された場合に1、そうでない場合に0となるバイナリ変数です。

　　制約条件

各点から1度だけ出る

\displaystyle \sum_{j=0, j \neq i}^{n} x_{ij} = 1 \quad \forall i \in \{0, 1, \ldots, n\}

各地点 $i$ からは、他のいずれかの地点 $j$ へ1度だけ移動する必要があります。

各点に1度だけ入る

\displaystyle \sum_{j=0, j \neq i}^{n} x_{ji} = 1 \quad \forall i \in \{0, 1, \ldots, n\}

各地点 $i$ には、他のいずれかの地点 $j$ から1度だけ移動してくる必要があります。

サブツアー除去制約

\displaystyle u_i - u_j + n x_{ij} \leq n - 1 \quad \forall i, j \in \{1, 2, \ldots, n\}, i \neq j

この制約は、サブツアー（部分的な巡回路）が発生しないようにするためのものです。 $u_i$ は補助変数で、地点 $i$ の順序を表します。

　Pythonによる実装

TSPを解くために、PuLPライブラリを使用します。

　　必要なモジュールの読み込み

まず、必要なモジュールを読み込みます。

以下、コードです。

import pulp

import numpy as np

import matplotlib.pyplot as plt

import networkx as nx

import pulp import numpy as np import matplotlib.pyplot as plt import networkx as nx

import pulp
import numpy as np
import matplotlib.pyplot as plt
import networkx as nx

　　サンプルデータの生成

サンプルデータを生成します。

以下、コードです。

# サンプルデータの生成

n = 5 # 顧客数

np.random.seed(42)

points = np.random.rand(n+1, 2) * 10 # 座標

# 距離行列の計算

dist = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(axis=-1))

# サンプルデータの生成 n = 5 # 顧客数 np.random.seed(42) points = np.random.rand(n+1, 2) * 10 # 座標 # 距離行列の計算 dist = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(axis=-1))

# サンプルデータの生成
n = 5  # 顧客数
np.random.seed(42)
points = np.random.rand(n+1, 2) * 10  # 座標

# 距離行列の計算
dist = np.sqrt(((points[:, None, :] - points[None, :, :]) ** 2).sum(axis=-1))

どのようなデータなのか見てみます。

以下、コードです。

# サンプルデータの表示

## ラベルの設定

labels = ['D'] + [str(i) for i in range(1, n+1)]

print('座標')

for i, point in enumerate(points):

print(f'{labels[i]}: {point}')

print('距離行列\n', dist)

# サンプルデータの表示 ## ラベルの設定 labels = ['D'] + [str(i) for i in range(1, n+1)] print('座標') for i, point in enumerate(points): print(f'{labels[i]}: {point}') print('距離行列\n', dist)

# サンプルデータの表示
## ラベルの設定
labels = ['D'] + [str(i) for i in range(1, n+1)]  
print('座標')
for i, point in enumerate(points):
    print(f'{labels[i]}: {point}')

print('距離行列\n', dist)

以下、実行結果です。

座標
D: [3.74540119 9.50714306]
1: [7.31993942 5.98658484]
2: [1.5601864 1.5599452]
3: [0.58083612 8.66176146]
4: [6.01115012 7.08072578]
5: [0.20584494 9.69909852]
距離行列
 [[0.         5.01713601 8.24215491 3.27553692 3.31980708 3.54475744]
 [5.01713601 0.         7.2642889  7.25066088 1.70589385 8.02453102]
 [8.24215491 7.2642889  0.         7.16902511 7.09155104 8.25106402]
 [3.27553692 7.25066088 7.16902511 0.         5.65579207 1.10303516]
 [3.31980708 1.70589385 7.09155104 5.65579207 0.         6.36847266]
 [3.54475744 8.02453102 8.25106402 1.10303516 6.36847266 0.        ]]

グラフで表現してみます。

以下、コードです。

import networkx as nx

# ネットワークのグラフ化

G = nx.DiGraph()

for i in range(n+1):

for j in range(n+1):

if i != j:

G.add_edge(i, j, weight=dist[i, j])

plt.figure(figsize=(12, 8))

pos = {i: points[i] for i in range(n+1)}

labels = {i: f"D" if i == 0 else str(i) for i in range(n+1)}

weights = nx.get_edge_attributes(G, 'weight')

nx.draw(

pos,

with_labels=True,

labels=labels,

node_size=1000,

node_color='red',

font_size=16,

font_color='white',

edge_color='blue')

nx.draw_networkx_edge_labels(

pos,

font_size=14,

edge_labels={

(i, j): f"{w:.2f}" for (i, j),

w in weights.items()})

plt.title("TSP problem graph")

plt.show()

import networkx as nx # ネットワークのグラフ化 G = nx.DiGraph() for i in range(n+1): for j in range(n+1): if i != j: G.add_edge(i, j, weight=dist[i, j]) plt.figure(figsize=(12, 8)) pos = {i: points[i] for i in range(n+1)} labels = {i: f"D" if i == 0 else str(i) for i in range(n+1)} weights = nx.get_edge_attributes(G, 'weight') nx.draw( G, pos, with_labels=True, labels=labels, node_size=1000, node_color='red', font_size=16, font_color='white', edge_color='blue') nx.draw_networkx_edge_labels( G, pos, font_size=14, edge_labels={ (i, j): f"{w:.2f}" for (i, j), w in weights.items()}) plt.title("TSP problem graph") plt.show()

import networkx as nx

# ネットワークのグラフ化
G = nx.DiGraph()
for i in range(n+1):
    for j in range(n+1):
        if i != j:
            G.add_edge(i, j, weight=dist[i, j])

plt.figure(figsize=(12, 8))
pos = {i: points[i] for i in range(n+1)}
labels = {i: f"D" if i == 0 else str(i) for i in range(n+1)}
weights = nx.get_edge_attributes(G, 'weight')
nx.draw(
    G, 
    pos, 
    with_labels=True, 
    labels=labels, 
    node_size=1000, 
    node_color='red', 
    font_size=16, 
    font_color='white', 
    edge_color='blue')
nx.draw_networkx_edge_labels(
    G, 
    pos, 
    font_size=14, 
    edge_labels={
        (i, j): f"{w:.2f}" for (i, j), 
        w in weights.items()})
plt.title("TSP problem graph")
plt.show()

以下、実行結果です。

　　最適化

最適化問題をモデル化し最適解を求めます。

以下、コードです。

# モデルの作成

prob = pulp.LpProblem("TSP", pulp.LpMinimize)

# 決定変数の定義

x = pulp.LpVariable.dicts(

"route",

((i, j) for i in range(n+1) for j in range(n+1) if i != j),

cat='Binary')

# 目的関数

prob += pulp.lpSum(

x[i, j] * dist[i, j] for i in range(n+1) for j in range(n+1) if i != j

)

# 制約条件

for i in range(n+1):

prob += pulp.lpSum(x[i, j] for j in range(n+1) if i != j) == 1 # 各点から1度だけ出る

prob += pulp.lpSum(x[j, i] for j in range(n+1) if i != j) == 1 # 各点に1度だけ入る

# サブツアー除去制約

u = pulp.LpVariable.dicts("u", range(1, n+1), lowBound=0, upBound=n)

for i in range(1, n+1):

for j in range(1, n+1):

if i != j:

prob += u[i] - u[j] + n * x[i, j] <= n - 1

# 問題を解く

prob.solve(pulp.PULP_CBC_CMD(msg=False))

# 結果の表示

print(f"最適な総移動距離: {pulp.value(prob.objective):.2f}")

# モデルの作成 prob = pulp.LpProblem("TSP", pulp.LpMinimize) # 決定変数の定義 x = pulp.LpVariable.dicts( "route", ((i, j) for i in range(n+1) for j in range(n+1) if i != j), cat='Binary') # 目的関数 prob += pulp.lpSum( x[i, j] * dist[i, j] for i in range(n+1) for j in range(n+1) if i != j ) # 制約条件 for i in range(n+1): prob += pulp.lpSum(x[i, j] for j in range(n+1) if i != j) == 1 # 各点から1度だけ出る prob += pulp.lpSum(x[j, i] for j in range(n+1) if i != j) == 1 # 各点に1度だけ入る # サブツアー除去制約 u = pulp.LpVariable.dicts("u", range(1, n+1), lowBound=0, upBound=n) for i in range(1, n+1): for j in range(1, n+1): if i != j: prob += u[i] - u[j] + n * x[i, j] <= n - 1 # 問題を解く prob.solve(pulp.PULP_CBC_CMD(msg=False)) # 結果の表示 print(f"最適な総移動距離: {pulp.value(prob.objective):.2f}")

# モデルの作成
prob = pulp.LpProblem("TSP", pulp.LpMinimize)

# 決定変数の定義
x = pulp.LpVariable.dicts(
    "route", 
    ((i, j) for i in range(n+1) for j in range(n+1) if i != j),
    cat='Binary')

# 目的関数
prob += pulp.lpSum(
    x[i, j] * dist[i, j] for i in range(n+1) for j in range(n+1) if i != j
)

# 制約条件
for i in range(n+1):
    prob += pulp.lpSum(x[i, j] for j in range(n+1) if i != j) == 1  # 各点から1度だけ出る
    prob += pulp.lpSum(x[j, i] for j in range(n+1) if i != j) == 1  # 各点に1度だけ入る

# サブツアー除去制約
u = pulp.LpVariable.dicts("u", range(1, n+1), lowBound=0, upBound=n)
for i in range(1, n+1):
    for j in range(1, n+1):
        if i != j:
            prob += u[i] - u[j] + n * x[i, j] <= n - 1

# 問題を解く
prob.solve(pulp.PULP_CBC_CMD(msg=False))

# 結果の表示
print(f"最適な総移動距離: {pulp.value(prob.objective):.2f}")

以下、実行結果です。

最適な総移動距離: 24.11

最適化の結果から、最適ルートを取得します。

以下、コードです。

# 最適ルートの取得

route = [0]

while len(route) < n + 2:

added = False

for j in range(n+1):

if j not in route and pulp.value(x[route[-1], j]) == 1:

route.append(j)

added = True

break

if not added:

break

print("最適ルート:", route)

# 最適ルートの取得 route = [0] while len(route) < n + 2: added = False for j in range(n+1): if j not in route and pulp.value(x[route[-1], j]) == 1: route.append(j) added = True break if not added: break print("最適ルート:", route)

# 最適ルートの取得
route = [0]
while len(route) < n + 2:
    added = False
    for j in range(n+1):
        if j not in route and pulp.value(x[route[-1], j]) == 1:
            route.append(j)
            added = True
            break
    if not added:
        break
        
print("最適ルート:", route)

以下、実行結果です。

最適ルート: [0, 4, 1, 2, 3, 5]

グラフで表現してみます。

以下、コードです。

# 結果の可視化

plt.figure(figsize=(12, 8))

plt.scatter(points[:, 0], points[:, 1], c='red', s=500)

plt.plot(points[route + [0], 0], points[route + [0], 1], c='blue')

for i, point in enumerate(points):

plt.annotate(

f"{'D' if i == 0 else i}",

(point[0], point[1]),

color='white',

fontsize=12,

ha='center', # フォントを真ん中寄せにする

va='center' # 垂直方向も真ん中寄せにする

)

plt.title("TSP optimal route", fontsize=12)

plt.xlabel("X", fontsize=15)

plt.ylabel("Y", fontsize=15)

plt.grid(True)

plt.show()

# 結果の可視化 plt.figure(figsize=(12, 8)) plt.scatter(points[:, 0], points[:, 1], c='red', s=500) plt.plot(points[route + [0], 0], points[route + [0], 1], c='blue') for i, point in enumerate(points): plt.annotate( f"{'D' if i == 0 else i}", (point[0], point[1]), color='white', fontsize=12, ha='center', # フォントを真ん中寄せにする va='center' # 垂直方向も真ん中寄せにする ) plt.title("TSP optimal route", fontsize=12) plt.xlabel("X", fontsize=15) plt.ylabel("Y", fontsize=15) plt.grid(True) plt.show()

# 結果の可視化
plt.figure(figsize=(12, 8))
plt.scatter(points[:, 0], points[:, 1], c='red', s=500)
plt.plot(points[route + [0], 0], points[route + [0], 1], c='blue')
for i, point in enumerate(points):
    plt.annotate(
        f"{'D' if i == 0 else i}",
        (point[0], point[1]), 
        color='white', 
        fontsize=12,
        ha='center',  # フォントを真ん中寄せにする
        va='center'   # 垂直方向も真ん中寄せにする
    )
plt.title("TSP optimal route", fontsize=12)
plt.xlabel("X", fontsize=15)
plt.ylabel("Y", fontsize=15)
plt.grid(True)
plt.show()

以下、実行結果です。

これらのコードを実行すると、以下のことが分かります。

最適ルート: 各顧客をどの順序で訪問するべきかが分かります。
総移動距離: このルートで移動した場合の総距離が計算されます。

この結果から、以下のような洞察が得られます。

最も効率的な顧客訪問順序
総移動距離に基づく燃料コストや所要時間の見積もり

TSPは配送ルート最適化の基本ですが、現実の配送問題はより複雑です。そこで登場するのが、車両経路問題（VRP: Vehicle Routing Problem）です。

車両経路問題（VRP）

　問題設定

VRPは、TSPを拡張し、複数の車両や容量制約、時間枠制約などを考慮した問題です。

具体的には、どの車両に、どの顧客を割り当て、どのような順番で顧客を回ればいいのか、その最短ルートを求める問題です。

要は、顧客の割り当て問題と、最短ルート問題を解く問題です。

与えられた情報

配送車両：K台
各車両の容量
顧客数：n
各顧客の需要量
各顧客に時間枠制約（要望時間）
各顧客のサービス時間

目的

総移動距離を最小化しつつ、全ての制約を満たす配送計画を立てる

決定変数

$x_{ijk}$ : 車両kが地点iから地点jへ移動する場合は1、そうでない場合は0

制約条件

各顧客（地点）を1回だけ訪問すること
各車両はデポから出発し、デポに戻ってくること
車両の容量制約を守ること
顧客（地点）の時間枠制約を守ること
時間の連続性（到着時間が出発時間や移動時間、サービス時間などを考慮）を守ること
部分巡回（サブツアー）の除去

　数理モデル

　　モデル全体

車両経路問題（VRP）を定式化すると、以下のようになります。

\begin{array}{ll}\text{Minimize} \quad & \displaystyle \sum_{k=1}^{K} \sum_{i=0}^{n} \sum_{j=0, j \neq i}^{n} d_{ij} x_{ijk} \\ \text{subject to} \quad & \\ & \displaystyle \sum_{i=0}^{n} x_{i0k} = 1 \quad \forall k \in \{1, \ldots, K\} \\ & \displaystyle \sum_{j=0}^{n} x_{0jk} = 1 \quad \forall k \in \{1, \ldots, K\} \\ & \displaystyle \sum_{k=1}^{K} \sum_{j=0, j \neq i}^{n} x_{ijk} = 1 \quad \forall i \in \{0, 1, \ldots, n\} \\ & \displaystyle \sum_{k=1}^{K} \sum_{i=0, i \neq j}^{n} x_{ijk} = 1 \quad \forall j \in \{0, 1, \ldots, n\} \\ & \displaystyle \sum_{k=1}^{K} \sum_{i=0}^{n} q_i x_{ijk} \leq Q_k \quad \forall k \in \{1, \ldots, K\} \\[0.5cm] & \displaystyle a_i \leq t_{ik} \leq b_i \quad \forall k \in \{1, \ldots, K\}, \forall i \in \{0, 1, \ldots, n\} \\[0.3cm] & \displaystyle t_{jk} \geq t_{ik} + s_i + d_{ij} - M (1 - x_{ijk}) \quad \forall k \in \{1, \ldots, K\}, \forall i, j \in \{0, 1, \ldots, n\}, i \neq j \\[0.3cm] & \displaystyle x_{ijk} \in \{0, 1\} \quad \forall i, j \in \{0, 1, \ldots, n\}, \forall k \in \{1, \ldots, K\} \\[0.3cm] & \displaystyle t_{ik} \geq 0 \quad \forall i \in \{0, 1, \ldots, n\}, \forall k \in \{1, \ldots, K\}\\[0.3cm] & \displaystyle u_i - u_j + n x_{ijk} \leq n - 1 \quad \forall k \in \{1, \ldots, K\}, \forall i, j \in \{1, \ldots, n\}, i \neq j\end{array}

以下、変数です。

$x_{ijk} \in \{0, 1\}$ : 車両 $k$ が地点 $i$ から地点 $j$ への移動を選択した場合に1、そうでない場合に0。
$t_{ik} \geq 0$ : 車両 $k$ が地点 $i$ に到着する時間。
$u_i \in \{0, 1, \ldots, n\}$ : 地点 $i$ の順序を表す補助変数。

　　目的関数

総移動距離を最小化することを目的としています。

\displaystyle \text{Minimize} \quad \sum_{k=1}^{K} \sum_{i=0}^{n} \sum_{j=0, j \neq i}^{n} d_{ij} x_{ijk}

ここで、 $d_{ij}$ は地点 $i$ から地点 $j$ への距離を表し、 $x_{ijk}$ は車両 $k$ が地点 $i$ から地点 $j$ への移動を選択した場合に1、そうでない場合に0となるバイナリ変数です。

　　制約条件

各車両はデポへ戻る

\displaystyle \sum_{i=0}^{n} x_{i0k} = 1 \quad \forall k \in \{1, \ldots, K\}

各車両 $k$ は、デポに帰ってきます。ただし、1度だけです。

各車両はデポから出る

\displaystyle \sum_{j=0}^{n} x_{0jk} = 1 \quad \forall k \in \{1, \ldots, K\}

各車両 $k$ は、デポから出発します。ただし、1度だけです。

いずれかの車両は各点（顧客）から1度だけ出る

\displaystyle \sum_{k=1}^{K} \sum_{j=0, j \neq i}^{n} x_{ijk} = 1 \quad \forall i \in \{0, 1, \ldots, n\}

各地点 $i$ からは、他のいずれかの地点 $j$ へ1度だけ移動する必要があります。

いずれかの車両は各点（顧客）に1度だけ入る

\displaystyle \sum_{k=1}^{K} \sum_{i=0, i \neq j}^{n} x_{ijk} = 1 \quad \forall j \in \{0, 1, \ldots, n\}

各地点 $j$ には、他のいずれかの地点 $i$ から1度だけ移動してくる必要があります。

車両の容量制約

\displaystyle \sum_{k=1}^{K} \sum_{i=0}^{n} q_i x_{ijk} \leq Q_k \quad \forall k \in \{1, \ldots, K\}

ここで、 $q_i$ は地点 $i$ の需要量を表し、 $Q_k$ は車両 $k$ の最大容量を表します。

時間枠制約

\displaystyle a_i \leq t_{ik} \leq b_i \quad \forall k \in \{1, \ldots, K\}, \forall i \in \{0, 1, \ldots, n\}

ここで、 $a_i$ と $b_i$ は地点 $i$ の時間枠の下限と上限を表し、 $t_{ik}$ は車両 $k$ が地点 $i$ に到着する時間を表します。

時間の連続性制約

\displaystyle t_{jk} \geq t_{ik} + s_i + d_{ij} - M (1 - x_{ijk}) \quad \forall k \in \{1, \ldots, K\}, \forall i, j \in \{0, 1, \ldots, n\}, i \neq j

ここで、 $s_i$ は地点 $i$ でのサービス時間を表し、 $M$ は十分に大きな定数です。この制約は、車両が次の地点に移動するための時間を考慮しています。

サブツアー除去制約

\displaystyle u_i - u_j + n x_{ijk} \leq n - 1 \quad \forall k \in \{1, \ldots, K\}, \forall i, j \in \{1, \ldots, n\}, i \neq j

この制約は、サブツアー（部分的な巡回路）が発生しないようにするためのものです。 $u_i$ は補助変数で、地点 $i$ の順序を表します。

　Pythonによる実装

VRPを解くために、OR-Toolsライブラリを使用します。OR-Toolsは、Googleが開発した最適化ソルバーで、VRPのような複雑な問題を効率的に解くことができます。

　　必要なモジュールの読み込み

まず、必要なモジュールを読み込みます。

以下、コードです。

# 必要なライブラリのインポート

from ortools.constraint_solver import routing_enums_pb2

from ortools.constraint_solver import pywrapcp

import numpy as np

import matplotlib.pyplot as plt

# 必要なライブラリのインポート from ortools.constraint_solver import routing_enums_pb2 from ortools.constraint_solver import pywrapcp import numpy as np import matplotlib.pyplot as plt

# 必要なライブラリのインポート
from ortools.constraint_solver import routing_enums_pb2
from ortools.constraint_solver import pywrapcp
import numpy as np
import matplotlib.pyplot as plt

車両経路問題（VRP: Vehicle Routing Problem）を解くために以下をインポートします。

from ortools.constraint_solver import routing_enums_pb2

ルーティングアルゴリズムの動作を制御するための定数を定義するためのものです。
具体的には、初期解法戦略や局所探索のメタヒューリスティックなどのアルゴリズムを設定するのに利用します。
初期解法戦略とは、ルーティング問題の初期解を生成するための方法のことです。
局所探索のメタヒューリスティックとは、初期解を改善するための探索アルゴリズムのことです。

from ortools.constraint_solver import pywrapcp

ルーティングモデルやインデックスマネージャーなど、ルーティング問題を解くためのクラスや関数です。
ルーティングモデルは、ルーティング問題の定式化と解法を管理するためのクラスです。
インデックスマネージャーは、ルーティングモデル内で使用されるノードのインデックスを管理するためのクラスです。
これらのコンポーネントを組み合わせ、ルーティング問題を定式化し解きます。

インデックスマネージャーとルーティングモデルのインスタンスを生成し、VRPの定式化などの設定を行います。

それをどのように解くのかを、初期解法戦略や局所探索のメタヒューリスティックなどのアルゴリズムを設定し、問題を解いていきます。

　　サンプルデータの生成

次に、サンプルデータを生成します。

このサンプルデータのデータモデルの構成要素は以下の通りです。

num_vehicles: 車両数
depot: デポ（出発地点）のインデックス
locations: 各地点の座標（デポを含む）
demands: 各地点の需要量（デポの需要は0）
vehicle_capacities: 各車両の容量
service_times: 各顧客のサービス時間
time_windows: 各地点の時間枠
distance_matrix: 各地点間の距離行列

では、生成します。サンプルデータは辞書型（dictionary）で格納します。

以下、コードです。

# データモデルの作成

def create_data_model():

"""問題のデータを格納します。"""

data = {}

data['num_vehicles'] = 3 # 車両数

data['depot'] = 0

# 10顧客 + デポのためのランダムな座標を生成

np.random.seed(42)

data['locations'] = np.random.rand(11, 2) * 100

# ランダムな需要を生成

data['demands'] = np.random.randint(1, 4, 11)

data['demands'][0] = 0 # デポの需要は0

# ランダムな車両容量を生成

data['vehicle_capacities'] = np.random.randint(10, 20, data['num_vehicles'])

# 顧客のサービス時間を設定（例: 一律10分）

data['service_times'] = np.full(11, 10)

data['service_times'][0] = 0 # デポのサービス時間は0

# 時間枠を生成

data['time_windows'] = [(0, 2000)] + [(np.random.randint(0, 200), np.random.randint(801, 2000)) for _ in range(10)]

# 距離行列を計算

data['distance_matrix'] = (np.sqrt(((data['locations'][:, None, :] - data['locations'][None, :, :]) ** 2).sum(axis=-1)) * 10).astype(int)

return data

# データモデルを作成

data = create_data_model()

# 確認用の出力

print("Vehicle capacities:", data['vehicle_capacities'])

print("Customer demands:", data['demands'])

print("Service times:", data['service_times'])

print("Time windows:", data['time_windows'])

# データモデルの作成 def create_data_model(): """問題のデータを格納します。""" data = {} data['num_vehicles'] = 3 # 車両数 data['depot'] = 0 # 10顧客 + デポのためのランダムな座標を生成 np.random.seed(42) data['locations'] = np.random.rand(11, 2) * 100 # ランダムな需要を生成 data['demands'] = np.random.randint(1, 4, 11) data['demands'][0] = 0 # デポの需要は0 # ランダムな車両容量を生成 data['vehicle_capacities'] = np.random.randint(10, 20, data['num_vehicles']) # 顧客のサービス時間を設定（例: 一律10分） data['service_times'] = np.full(11, 10) data['service_times'][0] = 0 # デポのサービス時間は0 # 時間枠を生成 data['time_windows'] = [(0, 2000)] + [(np.random.randint(0, 200), np.random.randint(801, 2000)) for _ in range(10)] # 距離行列を計算 data['distance_matrix'] = (np.sqrt(((data['locations'][:, None, :] - data['locations'][None, :, :]) ** 2).sum(axis=-1)) * 10).astype(int) return data # データモデルを作成 data = create_data_model() # 確認用の出力 print("Vehicle capacities:", data['vehicle_capacities']) print("Customer demands:", data['demands']) print("Service times:", data['service_times']) print("Time windows:", data['time_windows'])

# データモデルの作成
def create_data_model():
    """問題のデータを格納します。"""
    data = {}
    data['num_vehicles'] = 3  # 車両数
    data['depot'] = 0
    # 10顧客 + デポのためのランダムな座標を生成
    np.random.seed(42)
    data['locations'] = np.random.rand(11, 2) * 100
    # ランダムな需要を生成
    data['demands'] = np.random.randint(1, 4, 11) 
    data['demands'][0] = 0  # デポの需要は0
    # ランダムな車両容量を生成
    data['vehicle_capacities'] = np.random.randint(10, 20, data['num_vehicles'])
    # 顧客のサービス時間を設定（例: 一律10分）
    data['service_times'] = np.full(11, 10)
    data['service_times'][0] = 0  # デポのサービス時間は0
    # 時間枠を生成
    data['time_windows'] = [(0, 2000)] + [(np.random.randint(0, 200), np.random.randint(801, 2000)) for _ in range(10)]
    # 距離行列を計算
    data['distance_matrix'] = (np.sqrt(((data['locations'][:, None, :] - data['locations'][None, :, :]) ** 2).sum(axis=-1)) * 10).astype(int)
    return data

# データモデルを作成
data = create_data_model()

# 確認用の出力
print("Vehicle capacities:", data['vehicle_capacities'])
print("Customer demands:", data['demands'])
print("Service times:", data['service_times'])
print("Time windows:", data['time_windows'])

以下、実行結果です。

Vehicle capacities: [16 14 18]
Customer demands: [0 2 2 3 2 3 3 1 3 1 3]
Service times: [ 0 10 10 10 10 10 10 10 10 10 10]
Time windows: [(0, 2000), (166, 1188), (88, 1116), (13, 1042), (8, 1365), (129, 892), (110, 1756), (198, 1309), (7, 835), (80, 1826), (133, 1366)]

グラフで表現してみます。

以下、コードです。

import networkx as nx

# データモデルの取得

locations = data['locations']

distance_matrix = data['distance_matrix']

# グラフの作成

G = nx.DiGraph()

for i in range(len(locations)):

for j in range(len(locations)):

if i != j:

G.add_edge(i, j, weight=distance_matrix[i, j])

# グラフの描画

plt.figure(figsize=(12, 8))

pos = {i: locations[i] for i in range(len(locations))}

labels = {i: f"D" if i == 0 else str(i) for i in range(len(locations))}

weights = nx.get_edge_attributes(G, 'weight')

nx.draw(

pos,

with_labels=True,

labels=labels,

node_size=1000,

node_color='red',

font_size=16,

font_color='white',

edge_color='blue')

nx.draw_networkx_edge_labels(

pos,

font_size=14,

edge_labels={

(i, j): f"{w:.2f}" for (i, j),

w in weights.items()})

plt.title("VRP problem graph")

plt.show()

import networkx as nx # データモデルの取得 locations = data['locations'] distance_matrix = data['distance_matrix'] # グラフの作成 G = nx.DiGraph() for i in range(len(locations)): for j in range(len(locations)): if i != j: G.add_edge(i, j, weight=distance_matrix[i, j]) # グラフの描画 plt.figure(figsize=(12, 8)) pos = {i: locations[i] for i in range(len(locations))} labels = {i: f"D" if i == 0 else str(i) for i in range(len(locations))} weights = nx.get_edge_attributes(G, 'weight') nx.draw( G, pos, with_labels=True, labels=labels, node_size=1000, node_color='red', font_size=16, font_color='white', edge_color='blue') nx.draw_networkx_edge_labels( G, pos, font_size=14, edge_labels={ (i, j): f"{w:.2f}" for (i, j), w in weights.items()}) plt.title("VRP problem graph") plt.show()

import networkx as nx

# データモデルの取得
locations = data['locations']
distance_matrix = data['distance_matrix']

# グラフの作成
G = nx.DiGraph()
for i in range(len(locations)):
    for j in range(len(locations)):
        if i != j:
            G.add_edge(i, j, weight=distance_matrix[i, j])

# グラフの描画
plt.figure(figsize=(12, 8))
pos = {i: locations[i] for i in range(len(locations))}
labels = {i: f"D" if i == 0 else str(i) for i in range(len(locations))}
weights = nx.get_edge_attributes(G, 'weight')
nx.draw(
    G, 
    pos, 
    with_labels=True, 
    labels=labels, 
    node_size=1000, 
    node_color='red', 
    font_size=16, 
    font_color='white', 
    edge_color='blue')
nx.draw_networkx_edge_labels(
    G, 
    pos, 
    font_size=14, 
    edge_labels={
        (i, j): f"{w:.2f}" for (i, j), 
        w in weights.items()})
plt.title("VRP problem graph")
plt.show()

以下、実行結果です。

　　最適化

車両経路問題（VRP: Vehicle Routing Problem）を解くための準備を行います。

以下、コードです。

if solution:

plt.figure(figsize=(12, 8))

plt.scatter(data['locations'][:, 0], data['locations'][:, 1], c='red', s=500)

plt.scatter(data['locations'][0, 0], data['locations'][0, 1], c='green', s=100, marker='s')

for i, coord in enumerate(data['locations']):

plt.annotate(

str(i),

(coord[0], coord[1]),

fontsize=12,

color='white',

ha='center', # フォントを真ん中寄せにする

va='center' # 垂直方向も真ん中寄せにする

)

colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange']

for vehicle_id in range(data['num_vehicles']):

index = routing.Start(vehicle_id)

route = []

while not routing.IsEnd(index):

node_index = manager.IndexToNode(index)

route.append(node_index)

index = solution.Value(routing.NextVar(index))

route.append(manager.IndexToNode(index))

route_coords = data['locations'][route]

plt.plot(

route_coords[:, 0],

route_coords[:, 1],

c=colors[vehicle_id % len(colors)],

linewidth=2, label=f'Vehicle {vehicle_id}')

plt.title('VRP Solution')

plt.xlabel('X')

plt.ylabel('Y')

plt.legend()

plt.grid(True)

plt.show()

else:

print('No solution to visualize.')

if solution: plt.figure(figsize=(12, 8)) plt.scatter(data['locations'][:, 0], data['locations'][:, 1], c='red', s=500) plt.scatter(data['locations'][0, 0], data['locations'][0, 1], c='green', s=100, marker='s') for i, coord in enumerate(data['locations']): plt.annotate( str(i), (coord[0], coord[1]), fontsize=12, color='white', ha='center', # フォントを真ん中寄せにする va='center' # 垂直方向も真ん中寄せにする ) colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange'] for vehicle_id in range(data['num_vehicles']): index = routing.Start(vehicle_id) route = [] while not routing.IsEnd(index): node_index = manager.IndexToNode(index) route.append(node_index) index = solution.Value(routing.NextVar(index)) route.append(manager.IndexToNode(index)) route_coords = data['locations'][route] plt.plot( route_coords[:, 0], route_coords[:, 1], c=colors[vehicle_id % len(colors)], linewidth=2, label=f'Vehicle {vehicle_id}') plt.title('VRP Solution') plt.xlabel('X') plt.ylabel('Y') plt.legend() plt.grid(True) plt.show() else: print('No solution to visualize.')

if solution:
    plt.figure(figsize=(12, 8))
    plt.scatter(data['locations'][:, 0], data['locations'][:, 1], c='red', s=500)
    plt.scatter(data['locations'][0, 0], data['locations'][0, 1], c='green', s=100, marker='s') 
    for i, coord in enumerate(data['locations']):
        plt.annotate(
            str(i), 
            (coord[0], coord[1]), 
            fontsize=12, 
            color='white', 
            ha='center',  # フォントを真ん中寄せにする
            va='center'   # 垂直方向も真ん中寄せにする
        )
    
    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange']
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        route = []
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route.append(node_index)
            index = solution.Value(routing.NextVar(index))
        route.append(manager.IndexToNode(index))
        
        route_coords = data['locations'][route]
        plt.plot(
            route_coords[:, 0], 
            route_coords[:, 1], 
            c=colors[vehicle_id % len(colors)], 
            linewidth=2, label=f'Vehicle {vehicle_id}')
    
    plt.title('VRP Solution')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print('No solution to visualize.')

コードを簡単に説明します。

ルーティングインデックスマネージャーとモデルの作成

RoutingIndexManager
- ノード（地点）と車両のインデックスを管理するためのクラスです。
- これにより、ユーザーが定義したノードのインデックスと、ルーティングモデル内で使用される内部インデックスを相互に変換することができます。
RoutingModel
- ルーティング問題の定式化と解法を管理するためのクラスです。
- 具体的には、ノード間の移動コスト、制約条件（例えば、車両の容量や時間枠）、および目的関数（例えば、総移動距離の最小化）を定義します。

距離コールバックの追加

distance_callback関数は、2つのノード間の距離を計算して返す関数です。
この関数は、ルーティングモデルに登録され、各車両の移動コスト（距離）を計算するために使用されます。

容量制約の追加

demand_callback関数は、各ノードの需要量を計算して返す関数です。
この関数は、ルーティングモデルに登録され、各車両の容量制約を追加するために使用されます。
これにより、各車両がその容量を超えないように制約が追加されます。

時間枠制約の追加

time_callback関数は、2つのノード間の移動時間を計算して返す関数です。
この関数は、ルーティングモデルに登録され、各車両の時間枠制約を追加するために使用されます。
各ノードの時間枠を設定することで、特定の時間内に訪問する必要があるように制約が追加され、現実的なスケジュールを考慮したルートが生成されます。
これにより、車両経路問題を解くための準備が整いました。次のステップでは、探索パラメータを設定し、問題を解くことができます。

車両経路問題（VRP: Vehicle Routing Problem）を解きます。

以下、コードです。

# 初期解法の設定

search_parameters = pywrapcp.DefaultRoutingSearchParameters()

search_parameters.first_solution_strategy = (

routing_enums_pb2.FirstSolutionStrategy.AUTOMATIC)

search_parameters.local_search_metaheuristic = (

routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)

search_parameters.time_limit.FromSeconds(180)

# 問題を解く

solution = routing.SolveWithParameters(search_parameters)

# 初期解法の設定 search_parameters = pywrapcp.DefaultRoutingSearchParameters() search_parameters.first_solution_strategy = ( routing_enums_pb2.FirstSolutionStrategy.AUTOMATIC) search_parameters.local_search_metaheuristic = ( routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH) search_parameters.time_limit.FromSeconds(180) # 問題を解く solution = routing.SolveWithParameters(search_parameters)

# 初期解法の設定
search_parameters = pywrapcp.DefaultRoutingSearchParameters()
search_parameters.first_solution_strategy = (
    routing_enums_pb2.FirstSolutionStrategy.AUTOMATIC)
search_parameters.local_search_metaheuristic = (
    routing_enums_pb2.LocalSearchMetaheuristic.GUIDED_LOCAL_SEARCH)
search_parameters.time_limit.FromSeconds(180) 

# 問題を解く
solution = routing.SolveWithParameters(search_parameters)

このコードで、車両経路問題（VRP: Vehicle Routing Problem）を解くための探索パラメータを設定し、問題を解きます。

初期解法の設定

search_parametersのインスタンスを作成します。
first_solution_strategyをAUTOMATICに設定します。これは、初期解法の戦略を自動的に選択することを意味します。
local_search_metaheuristicをGUIDED_LOCAL_SEARCHに設定します。これは、局所探索のメタヒューリスティックとしてガイド付き局所探索アルゴリズムを使用することを意味します。
time_limitを180秒に設定します。これは、探索にかける最大時間を180秒に制限することを意味します。

問題を解く

routing.SolveWithParametersメソッドを使用して、設定した探索パラメータで車両経路問題を解きます。
解が見つかった場合、solutionに解が格納されます。

では、解を表示します。

以下、コードです。

# 解の表示

if solution:

total_distance = 0

total_load = 0

for vehicle_id in range(data['num_vehicles']):

index = routing.Start(vehicle_id)

plan_output = f'Route for vehicle {vehicle_id}:\n'

route_distance = 0

route_load = 0

while not routing.IsEnd(index):

node_index = manager.IndexToNode(index)

route_load += data['demands'][node_index]

plan_output += f' {node_index} Load({route_load}) -> '

previous_index = index

index = solution.Value(routing.NextVar(index))

route_distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id)

plan_output += f' {manager.IndexToNode(index)} Load({route_load})\n'

plan_output += f'Distance of the route: {route_distance}m\n'

plan_output += f'Load of the route: {route_load}\n'

print(plan_output)

total_distance += route_distance

total_load += route_load

print(f'Total distance of all routes: {total_distance}m')

print(f'Total load of all routes: {total_load}')

else:

print('No solution found !')

# 解の表示 if solution: total_distance = 0 total_load = 0 for vehicle_id in range(data['num_vehicles']): index = routing.Start(vehicle_id) plan_output = f'Route for vehicle {vehicle_id}:\n' route_distance = 0 route_load = 0 while not routing.IsEnd(index): node_index = manager.IndexToNode(index) route_load += data['demands'][node_index] plan_output += f' {node_index} Load({route_load}) -> ' previous_index = index index = solution.Value(routing.NextVar(index)) route_distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id) plan_output += f' {manager.IndexToNode(index)} Load({route_load})\n' plan_output += f'Distance of the route: {route_distance}m\n' plan_output += f'Load of the route: {route_load}\n' print(plan_output) total_distance += route_distance total_load += route_load print(f'Total distance of all routes: {total_distance}m') print(f'Total load of all routes: {total_load}') else: print('No solution found !')

# 解の表示
if solution:
    total_distance = 0
    total_load = 0
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        plan_output = f'Route for vehicle {vehicle_id}:\n'
        route_distance = 0
        route_load = 0
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route_load += data['demands'][node_index]
            plan_output += f' {node_index} Load({route_load}) -> '
            previous_index = index
            index = solution.Value(routing.NextVar(index))
            route_distance += routing.GetArcCostForVehicle(previous_index, index, vehicle_id)
        plan_output += f' {manager.IndexToNode(index)} Load({route_load})\n'
        plan_output += f'Distance of the route: {route_distance}m\n'
        plan_output += f'Load of the route: {route_load}\n'
        print(plan_output)
        total_distance += route_distance
        total_load += route_load
    print(f'Total distance of all routes: {total_distance}m')
    print(f'Total load of all routes: {total_load}')
else:
    print('No solution found !')

以下、実行結果です。

Route for vehicle 0:
 0 Load(0) ->  8 Load(3) ->  7 Load(4) ->  2 Load(6) ->  0 Load(6)
Distance of the route: 1654m
Load of the route: 6

Route for vehicle 1:
 0 Load(0) ->  4 Load(2) ->  1 Load(4) ->  6 Load(7) ->  10 Load(10) ->  9 Load(11) ->  0 Load(11)
Distance of the route: 2028m
Load of the route: 11

Route for vehicle 2:
 0 Load(0) ->  3 Load(3) ->  5 Load(6) ->  0 Load(6)
Distance of the route: 791m
Load of the route: 6

Total distance of all routes: 4473m
Total load of all routes: 23

グラフで表現してみます。

以下、コードです。

if solution:

plt.figure(figsize=(12, 8))

plt.scatter(data['locations'][:, 0], data['locations'][:, 1], c='red', s=500)

plt.scatter(data['locations'][0, 0], data['locations'][0, 1], c='green', s=100, marker='s')

for i, coord in enumerate(data['locations']):

plt.annotate(

str(i),

(coord[0], coord[1]),

fontsize=12,

color='white',

ha='center', # フォントを真ん中寄せにする

va='center' # 垂直方向も真ん中寄せにする

)

colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange']

for vehicle_id in range(data['num_vehicles']):

index = routing.Start(vehicle_id)

route = []

while not routing.IsEnd(index):

node_index = manager.IndexToNode(index)

route.append(node_index)

index = solution.Value(routing.NextVar(index))

route.append(manager.IndexToNode(index))

route_coords = data['locations'][route]

plt.plot(

route_coords[:, 0],

route_coords[:, 1],

c=colors[vehicle_id % len(colors)],

linewidth=2, label=f'Vehicle {vehicle_id}')

plt.title('VRP Solution')

plt.xlabel('X')

plt.ylabel('Y')

plt.legend()

plt.grid(True)

plt.show()

else:

print('No solution to visualize.')

if solution:
    plt.figure(figsize=(12, 8))
    plt.scatter(data['locations'][:, 0], data['locations'][:, 1], c='red', s=500)
    plt.scatter(data['locations'][0, 0], data['locations'][0, 1], c='green', s=100, marker='s') 
    for i, coord in enumerate(data['locations']):
        plt.annotate(
            str(i), 
            (coord[0], coord[1]), 
            fontsize=12, 
            color='white', 
            ha='center',  # フォントを真ん中寄せにする
            va='center'   # 垂直方向も真ん中寄せにする
        )
    
    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k', 'orange']
    for vehicle_id in range(data['num_vehicles']):
        index = routing.Start(vehicle_id)
        route = []
        while not routing.IsEnd(index):
            node_index = manager.IndexToNode(index)
            route.append(node_index)
            index = solution.Value(routing.NextVar(index))
        route.append(manager.IndexToNode(index))
        
        route_coords = data['locations'][route]
        plt.plot(
            route_coords[:, 0], 
            route_coords[:, 1], 
            c=colors[vehicle_id % len(colors)], 
            linewidth=2, label=f'Vehicle {vehicle_id}')
    
    plt.title('VRP Solution')
    plt.xlabel('X')
    plt.ylabel('Y')
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print('No solution to visualize.')

以下、実行結果です。