容量制約つき運搬経路問題(Capacitated Vehicle Routing Problem, CVRP)¶
このチュートリアルでは、運搬経路問題 (VRP) の一種である容量制約つき運搬経路問題を扱います。
具体的な応用先として、
- 郵便などの運送業における効率的な配送計画の策定
- ごみ収集や道路清掃における訪問順序の決定
などがあります。
運搬経路問題とは、配送拠点 (depot) から複数の都市への配送を効率的に行おうとする配送ルート決定問題です。より具体的には、配送車両の総移動距離や総コストなどが最小になるような配送車両と都市の割り当て、都市の訪問順序を決定します。運搬経路問題は巡回セールスマン問題の一般化だと解釈することもできます。
本チュートリアルで取り扱う容量制約付き運搬経路問題は、上記運搬経路問題に各車両の積載上限が追加された問題です。つまり、各配送車両は積載量制約を満たした上で配送を行う必要があります。この問題は「巡回セールスマン問題+ナップサック問題」のように解釈することもできます。
今回は配送拠点(デポ)が一つかつ、各都市の需要と車の容量が整数値のみを取るような場合を考えます。
定式化¶
まずは定式化に必要な定数・変数を定義します。
定数¶
- $N$:都市の数
- $L$:積載可能量を超えない範囲で巡回可能な都市の最大数+2(出発点+到着点)
- $K$:車両の数
- $Q$:各車両の積載可能量、整数値
- $w_i$:都市 $i$ に配送する荷物の重さを表す整数
- $D_{i\to j}$:都市 $i$ から都市 $j$ までの距離を表現する実数(都市 $0$ = デポ)($i, j \in \{0,1,\dots,N \}$)
変数¶
- $x_{i,j,k} \in \{0, 1\} \quad (i \in \{1, \dots , L\}, j \in \{0, \dots, N\}, k \in \{1, \dots,
K\})$
車両$k$が、$i$番目に訪れる都市として都市 $j$ を選択するかどうかを表すバイナリ変数($x_{i,j,k}=1\Leftrightarrow$ 車両 $k$ が、$i$ 番目に訪れる都市として都市 $j$ を選択する)
続いて、変数が満たすべき制約を考えます。
制約条件¶
必要な制約は以下の4つです。
- 全ての車両はデポ(都市 $0$)から出発し、デポに到着する
- 車両 $k$ が $i$ 番目に訪問する場所は一箇所だけ
- 都市 $j$(デポを除く)はいずれかの車両によってちょうど一回だけ訪問される
- 車両 $k$ が運ぶ荷物の総重量は $Q$ 以下
これらの条件を式にすると以下のようになります。
-
全ての車両はデポ(都市 $0$)から出発し、デポに到着する
$ \begin{align} x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \end{align} $ -
車両 $k$ が $i$ 番目に訪問する場所は一箇所だけ
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$ -
都市 $j$(デポを除く)はいずれかの車両によってちょうど一回だけ訪問される
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$ -
車両 $k$ が運ぶ荷物の総重量は $Q$ 以下
$$ \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\}) $$
目的関数¶
今回は車両の総移動距離を最小化することを目指します。車両の総移動距離を数式で表現すると以下のようになります。 $$ \sum_{k=1}^K\sum_{i=1}^L\sum_{j_1=0}^N\sum_{j_2=0}^N D_{j_1\to j_2}x_{i,j_1,k}x_{i+1,j_2,k} $$
具体例を使った制約式の説明¶
ここでは具体例を使ってより感覚的な説明を行います。
以下は、制約条件を満たす変数の割当ての一例を、車両ごとの表にまとめたものです。横軸はアルファベットの都市名(depot は出発地)、縦軸はその都市を何番目に訪問するか(0 は depot を出発、9 は depot に到着)を表します。表のマスが $1$ のとき、その車両は都市に訪れることを意味します。
- 車両1 ($k=1$)
$x_{i,j,1}$ | depot | A | B | C | D | E | F | G | H | I |
---|---|---|---|---|---|---|---|---|---|---|
0(出発) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 |
3 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
4 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
5 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
6 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
7 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
8 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
9(到着) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
- 車両2 ($k=2$)
$x_{i,j,2}$ | depot | A | B | C | D | E | F | G | H | I |
---|---|---|---|---|---|---|---|---|---|---|
0(出発) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
2 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
3 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
4 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 |
6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |
7 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
8 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
9(到着) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
- 車両3 ($k=3$)
$x_{i,j,3}$ | depot | A | B | C | D | E | F | G | H | I |
---|---|---|---|---|---|---|---|---|---|---|
0(出発) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
1 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 |
2 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 |
3 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
4 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
5 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
6 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
7 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
8 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
9(到着) | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
先の『制約条件2: 車両 $k$ が $i$ 番目に訪問する場所は一箇所だけ』は、各行(横のライン)に現れる $1$ が一つだけになることと対応しており、次式で表現することができます。
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$
また、『制約条件3: 都市 $j$ はいずれかの車両によってちょうど一回だけ訪問される』は、全ての車両を合わせた、depot を除く各列(縦のライン)に現れる $1$ が一つだけになることと対応しており、次式で表現することができます。
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$
これらの制約によって全ての都市がいずれかの車両によってちょうど一回だけ訪問されることが保証されますが、配送拠点 (depot) は何回でも訪問可能な定式化になっています。そのため、各車両が配送開始時と配送終了時の二回だけ配送拠点を訪問するように後処理を行います。具体的には、表中の0行目と9行目以外で depot の列に $1$ が立っている場合、それを無視する($0$ だとみなす)という操作です。この操作によって得られる運搬経路に対応する移動距離が操作前のものと同じ、もしくは短くなるということが三角不等式を用いて証明できます。
import numpy as np
from geopy.distance import geodesic
# 車両数
nvehicle = 2
# 都市の数
ncity = 10
dimension = ncity + 1
avg_cities_per_vehicle = dimension // nvehicle
# 各都市における配送需要(重量)を乱数で決定
demand = np.random.randint(1, 100, size=dimension)
demand[0] = 0
demand_max = max(demand)
demand_mean = sum(demand) / len(demand)
# 全体的な需要に合わせ、車両の積載可能量 Q を設定する。
capacity = demand_max + int(demand_mean) * avg_cities_per_vehicle
次にそれぞれの都市の座標を乱数を使って生成します。全ての都市 $i,j$ 間の座標距離を計算することによって $D_{i \to j}$ を生成します。
# 座標の取り得る範囲を設定
lat_range = [35.7014, 35.968]
lon_range = [139.34, 140.04]
# 各都市の座標
ind2coord = {
i: (
np.random.uniform(lon_range[0], lon_range[1]),
np.random.uniform(lat_range[0], lat_range[1]),
)
for i in range(dimension)
}
# 2都市間の座標距離行列 D
distance_matrix = np.array(
[
[geodesic(coord_i[::-1], coord_j[::-1]).m for coord_j in ind2coord.values()]
for coord_i in ind2coord.values()
]
)
問題のプロット時に使用する色のリストを作成します。
_colors = [
"green",
"orange",
"blue",
"red",
"purple",
"pink",
"darkblue",
"cadetblue",
"darkred",
"lightred",
"darkgreen",
"lightgreen",
"lightblue",
"darkpurple",
]
都市の座標をプロットする関数 plot_solution
を定義します。
import folium
def plot_solution(coord: dict, title: str, best_tour: dict = dict()):
l = len(coord)
center = [
sum(lat for _, lat in coord.values()) / l,
sum(lon for lon, _ in coord.values()) / l,
]
m = folium.Map(center, tiles="Stamen Toner", zoom_start=10.5)
folium.Marker(
location=coord[0][::-1],
popup=f"depot",
icon=folium.Icon(icon="car", prefix="fa"),
).add_to(m)
_color = _colors[1]
if best_tour:
for k, tour in best_tour.items():
_color = _colors[k % len(_colors)]
for city in tour:
if city == 0:
continue
folium.Marker(
location=coord[city][::-1],
popup=f"person{k}",
icon=folium.Icon(
icon="user", prefix="fa", color="white", icon_color=_color
),
).add_to(m)
folium.vector_layers.PolyLine(
locations=[coord[city][::-1] for city in tour], color=_color, weight=3
).add_to(m)
else:
for k, node in coord.items():
if k == 0:
continue
folium.Marker(
location=node[::-1],
popup=f"customer{k}",
icon=folium.Icon(
icon="user", prefix="fa", color="white", icon_color=_color
),
).add_to(m)
title = f"<h4>{title}</h4>"
m.get_root().html.add_child(folium.Element(title))
return m
作成した問題を plot_solution
関数で表示します。車型のピンは depot の位置(出発点&到着点)、人型のピンは通るべき地点を表します。
# 埼玉・千葉・東京近辺
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}"
plot_solution(ind2coord, title)
訪問都市の上限の計算¶
決定変数のサイズを削減するための工夫として、一台の車が訪問できる都市の数の上限をあらかじめ計算しておきます。
以下の upperbound_of_tour
関数は、積載可能量を超えない範囲で最大いくつの都市を巡回できるかを計算する関数です。具体的には次のステップ計算します。
- まだ選んでいない都市の中で最も配送需要(重量)が少ない都市を選び、最大積載量からその都市の需要を減算する
- 最大積載量が0未満になるまでステップを繰り返し、繰り返した回数を戻り値とする
def upperbound_of_tour(capacity: int, demand: np.ndarray) -> int:
max_tourable_cities = 1
for w in sorted(demand):
capacity -= w
if capacity >= 0:
max_tourable_cities += 1
else:
return max_tourable_cities
return max_tourable_cities
二次多項式モデルの構築¶
次に、必要な決定変数を定義します。Amplify の BinarySymbolGenerator
を用います。合計$N+1$ 個の都市及びデポに対して、最大 $L$
都市の訪問を、$K$ 台の車両で行うので、$L \times (N+1) \times K$ の三次元配列としてQUBO変数を次のように定義します。
from amplify import BinaryPoly, BinarySymbolGenerator
gen = BinarySymbolGenerator()
# 積載可能量から1台の車両が訪問できる都市の最大数
max_tourable_cities = upperbound_of_tour(capacity, demand)
x = gen.array(max_tourable_cities, dimension, nvehicle)
前述の通り、制約条件は、
- 全ての車両はデポ(都市 $0$)から出発し、デポに到着する
- 車両 $k$ が $i$ 番目に訪問する場所は一箇所だけ
- 都市 $j$(デポを除く)はいずれかの車両によってちょうど一回だけ訪問される
- 車両 $k$ が運ぶ荷物の総重量は $Q$ 以下
つまり、
$$ \begin{align*} (1) \quad & x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ & x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ (2)\quad & \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\})\\ (3)\quad & \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\})\\ (4)\quad & \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\})\\ \end{align*} $$ です。
今回は less_equal
関数を用いて不等式制約を表現し、 equal_to
関数で等式制約を表現します。また、制約項の強さを決定する係数は、$\max
d_{ij}$ として実装しています。これにより、少なくとも制約項より目的関数が優先される事態を避けることができます。
制約条件①: 全ての車両はデポ(都市 0)から出発し、デポに到着する¶
制約式は以下で表されます。
$$ \begin{align*} x_{1,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \\ x_{L,j,k} = \begin{cases} 1 \quad (j = 0) \\ 0 \quad (j \neq 0) \end{cases} \end{align*} $$
この制約については代入によって実現します。代入することで求解の対象となる変数を節約することが出来ます。次のように直接値を代入する関数を作成します。
変数の値を事前に固定するには、x[i][j][k] = 1
のようにします。決定変数への代入操作は目的関数や制約式を定義する前に行う必要があります。
x[0][1:][:] = 0
x[-1][1:][:] = 0
x[0][0][:] = 1
x[-1][0][:] = 1
制約条件②: 車両$k$が$i$番目に訪問する場所は一箇所だけ¶
本制約条件は、次式で表されます。
$$ \sum_{j=0}^N x_{i,j,k} = 1 \quad (i \in \{1, \dots, L\}, k \in \{1, \dots, K\}) $$
本制約式を次のように実装します。
from amplify.constraint import one_hot
from amplify import sum_poly
def make_onetrip_constraints(x) -> list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
one_hot(sum_poly(dimension, lambda j: x[i][j][k]))
for i in range(max_tourable_cities)
for k in range(nvehicle)
]
return constraints
制約条件③: 都市$j$はいずれかの車両によってちょうど一回だけ訪問される¶
本制約条件は次式で表されます。
$$ \sum_{i=1}^L \sum_{k=1}^K x_{i,j,k}=1 \quad (j \in \{1, \dots, N\}) $$
これを実装すると、次のようになります。
def make_onevisit_constraints(x) -> list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
one_hot(
sum_poly(
max_tourable_cities, lambda i: sum_poly(nvehicle, lambda k: x[i][j][k])
)
)
for j in range(1, dimension)
]
return constraints
制約条件④: 車両 $k$ が運ぶ荷物の総重量は $Q$ 以下¶
本制約条件は次式で表されます。
$$ \sum_{i=1}^L \sum_{j=1}^N w_j x_{i,j,k} \leq Q \quad (k \in \{1, \dots, K\}) $$
この制約式は、less_equal
を使って実装することができます。また、今回の問題設定では必ずしも必要ではありませんが、less_equal
の
method
引数に InequalityFormulation.RelaxationQuadra
を指定します。制約式に対する厳密なペナルティ関数ではなくなりますが、多くの補助変数を節約するとともに、実数係数にも対応でき、問題サイズが大きい条件における求解可能性を上げることができます。
詳しくは https://amplify.fixstars.com/ja/docs/constraint.html#constraint-inequality-relaxation を参照してください。
from amplify.constraint import less_equal
from amplify import InequalityFormulation
def make_capacity_constraints(x, demand: np.ndarray, capacity: int) -> list:
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
constraints = [
less_equal(
demand * x[:, :, k],
capacity,
method=InequalityFormulation.RelaxationQuadra,
)
/ capacity
/ capacity
for k in range(nvehicle)
]
return constraints
from amplify import einsum
max_tourable_cities = x.shape[0]
dimension = x.shape[1]
nvehicle = x.shape[2]
xr = x.roll(-1, axis=0)
# 経路の総距離
objective = einsum("pq,ipk,iqk", distance_matrix, x, xr)
ここまでに実装した代入に基づく1つの制約条件、制約式に基づく3つの制約条件及び目的関数に基づき、次のように CVRP の QUBO
モデルを構築することができます。ここで、np.max(distance_matrix)
制約式に対する重みを設定しています。
constraints1 = make_onetrip_constraints(x)
constraints2 = make_onevisit_constraints(x)
constraints3 = make_capacity_constraints(x, demand, capacity)
constraints = sum(constraints1) + sum(constraints2) + sum(constraints3)
model = objective + constraints * np.max(distance_matrix)
Amplify クライアントの設定¶
Amplify のイジングマシン・量子アニーリングのクライアントを作成します。
from amplify.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 2000 # タイムアウト2秒
# client.token = "API トークンを入力してください"
作成したクライアントで Solver
をインスタンス化し、QUBOモデルに基づき求解を行います。
from amplify import Solver
solver = Solver(client)
result = solver.solve(model)
if len(result.solutions) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
x_values = result.solutions[0].values # 変数値の取り出し
def onehot2_sequence(x_values) -> dict[int, list]:
nvehicle = x_values.shape[2]
sequence = dict()
for k in range(nvehicle):
sequence[k] = np.where(x_values[:, :, k])[1]
return sequence
onehot2_sequence
関数で作ったリストからデポへの余計な訪問を取り除いたリストを返す process_sequence
関数を作成します。
def process_sequence(sequence: dict[int, list]) -> dict[int, list]:
new_seq = dict()
for k, v in sequence.items():
v = np.append(v, v[0])
mask = np.concatenate(([True], np.diff(v) != 0))
new_seq[k] = v[mask]
return new_seq
それでは、実際に結果を取り出し、可視化してみましょう。
変数配列 x
の decode
メンバ関数に解である x_values
を与え、変数配列に結果を代入します。
その後、先程定義した関数 onehot2_sequence
、process_sequence
を使って後処理を行います。
solution = x.decode(x_values)
sequence = onehot2_sequence(
solution
) # one-hot な変数テーブルを辞書に変換。key:車両インデックス, value:各車両が訪問した都市の順番が入ったリスト
best_tour = process_sequence(sequence) # 上の辞書からデポへの余計な訪問を取り除く
print(result.solutions[0].energy) # 目的関数値を表示
print(*best_tour.items(), sep="\n") # 求めた解を表示
plot_solution
関数で結果の可視化を行います。地図上にそれぞれの車両が通るべき経路が表示されます。
title = f"capacity={capacity}, ncity={ncity}, nvehicle={nvehicle}, cost={result.solutions[0].energy:.2f}"
plot_solution(ind2coord, title, best_tour)