配送計画問題
配送計画問題（運搬経路問題）は，我が国では最も普及しているロジスティクス・ツールである配送計画の基幹を成すモデルである． 本来ならば，配送だけでなく集荷にも使われるので運搬経路問題（vehicle routing problem）とよぶのが学術用語としては正しい使い方だが， 「運搬」という言葉のイメージが悪いためか，実務家および研究者の間でも配送計画問題とよばれることが多いので，ここでもそれにならうものとする．
一般に，配送計画問題の基本形は以下の仮定をもつ．
- デポとよばれる特定の地点を出発した運搬車が，顧客を経由し再びデポに戻る．このとき運搬車による顧客の通過順をルートとよぶ．
- デポに待機している運搬車の種類および最大積載重量は既知である．
- 顧客の位置は既知であり，各顧客の需要量も事前に与えられている．
- 地点間の移動時間，移動距離，移動費用は既知である．
- 1つのルートに含まれる顧客の需要量の合計は運搬車の最大積載重量を超えない．
- 運搬車の台数は，決められた上限を超えない．（超過した運搬車に対するレンタル料を考える場合もある．）
- 運搬車の稼働時間が与えられた上限を超えない．（超過時間を残業費用として考える場合もある．）
配送計画問題の応用としては，小売店への配送計画，スクールバスの巡回路決定， 郵便や宅配便の配達，ゴミの収集，燃料の配送などがある． もちろん，これらの応用に適用する際には，上の基本条件に新たな条件を付加する必要がある．
配送計画問題に対する厳密解法（最適解を算出することが保証されているアルゴリズム）としては列生成法（分枝価格法）が有効である． 列生成法は，時間枠制約（何時から何時までの間に顧客を訪問をしなければならないという制約条件） がきついときや，$1$台の運搬車が訪問する顧客数が少ないときに有効になる．
一方，一般の配送計画問題に対しては，現状ではヒューリスティクス（最適解を算出することが保証されていないアルゴリズム）を用いることが多い．
対称容量制約付き配送計画問題
運搬車の数を $m$，点（顧客およびデポを表す）の数を $n$ とする． デポの番号は $1$ とする． 顧客 $i=2,3,\ldots,n$ は需要量 $q_i$ をもち，その需要はある運搬車によって運ばれる（または収集される） ものとする．運搬車 $k=1,2,\ldots,m$ は有限の積載量上限 $Q$ をもち， 運搬車によって運ばれる需要量の合計は，その値を超えないものとする． 通常は，顧客の需要量の最大値 $\max \{ q_i \}$ は，運搬車の容量 $Q$ を超えないものと仮定する． もし，積載量の最大値を超える需要をもつ顧客が存在するなら， （積載量の上限に収まるように）需要を適当に分割して おくことによって，上の仮定を満たすように変換できる．
運搬車が点 $i$ から点 $j$ に移動するときにかかる費用を $c_{ij}$ と書く．ここでは移動費用は対称 (すなわち $c_{ij}=c_{ji}$) とする． 配送計画問題の目的はすべての顧客の需要を満たす $m$ 台の運搬車の 最適ルート（デポを出発して再びデポへ戻ってくる単純閉路）を求めることである．
点 $i,j$ 間を運搬車が移動する回数を表す変数 $x_{ij}$ を導入する． 対称な問題を仮定するので，変数 $x_{ij}$ は $i<j$ を満たす点 $i,j$ 間にだけ定義される． $x_{ij}$ がデポに接続しない枝に対しては，運搬車が通過するときに $1$，それ以外のとき $0$ を表すが， デポからある点 $j$ に移動してすぐにデポに帰還するいわゆるピストン輸送の場合には， $x_{1j}$ は $2$ となる．
容量制約付き配送計画問題の定式化は，以下のようになる．
$$ \begin{array}{l l l} minimize & \sum_{i,j} c_{ij} x_{ij} & \\ s.t. & \sum_{j} x_{1j} = 2 m & \\ & \sum_{j} x_{ij} = 2 & \forall i=2,3,\ldots,n \\ &\sum_{i,j \in S} x_{ij} \leq |S|-N(S) & \forall S \subset \{2,3,\ldots,n\}, |S| \geq 3 \\ & x_{1j} \in \{0,1,2\} & \forall j=2,\ldots,n \\ & x_{ij} \in \{0,1\} & \forall i <j: i \neq 1 \end{array} $$ここで，最初の制約は，デポ（点 $1$ ）から $m$台の運搬車があることを規定する． すなわち，点$1$ に出入りする運搬車を表す枝の本数が $2m$ 本であることを表す． 2番目の制約は，各顧客に1台の運搬車が訪問することを表す． 3番目の制約は，運搬車の容量制約と部分巡回路を禁止することを同時に規定する制約である． この制約で用いられる $N(S)$ は，顧客の部分集合 $S$ を与えたときに計算される関数であり， 「$S$ 内の顧客の需要を運ぶために必要な運搬車の数」と定義される．
$N(S)$ を計算するためには，箱詰め問題を解く必要があるが，通常は以下の下界で代用する． $$ \left\lceil \sum_{i \in S} q_i /Q \right\rceil $$
ここで $\lceil \cdot \rceil$ は天井関数である． 巡回セールスマン問題に対して適用したように， 線形緩和問題の解を $\bar{x}_e (e \in E)$ としたとき， $\bar{x}_e$ が正の枝から構成されるグラフ上で連結成分を求めることによる分枝カット法を考える．
容量制約付き配送計画問題のベンチマーク問題例は，以下のサイトからダウンロードできる．
https://neo.lcc.uma.es/vrp/vrp-instances/capacitated-vrp-instances/
ベンチマーク問題例を読み込む．
def distance(x1, y1, x2, y2):
"""distance: euclidean distance between (x1,y1) and (x2,y2)"""
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
folder = "../data/cvrp/"
# file_name = "E-n22-k4.vrp"
file_name = "E-n30-k3.vrp"
# file_name = "E-n51-k5.vrp"
f = open(folder + file_name)
data = f.readlines()
f.close()
n = int(data[3].split()[-1])
Q = int(data[5].split()[-1])
print("n=", n, "Q=", Q)
x, y, q = {}, {}, {}
pos = {}
for i, row in enumerate(data[7 : 7 + n]):
id_no, x[i], y[i] = list(map(int, row.split()))
pos[i] = x[i], y[i]
for i, row in enumerate(data[8 + n : 8 + 2 * n]):
id_no, q[i] = list(map(int, row.split()))
m = 3
c = {}
for i in range(n):
for j in range(n):
if i != j:
c[i, j] = int(distance(x[i], y[i], x[j], y[j]))
V = list(range(n))
model, vrp_callback = vrp(V, c, m, q, Q)
# model.Params.OutputFlag = 0 # silent mode
model.params.DualReductions = 0
model.params.LazyConstraints = 1
model.optimize(vrp_callback)
xx = model.__data
edges = []
for (i, j) in xx:
if xx[i, j].X > 0.5:
edges.append((i, j))
G = nx.Graph()
G.add_edges_from(edges)
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
非対称容量制約付き配送計画問題
非対称な移動費用をもつ容量制約付きの配送計画問題に対しては， 巡回セールスマン問題に対するMiller--Tucker--Zemlinのポテンシャル定式化を拡張したものが使える．
$$ \begin{array}{lll} minimize & \sum_{i,j} c_{ij} x_{ij} & \\ s.t. & \sum_{j} x_{ji} = \sum_{j} x_{ij} = 1 & \forall i \\ & \sum_{j} x_{0j} = \sum_{j} x_{j0} = m & \\ & u_i - u_j + Q x_{ij} \leq Q -q_j & \forall i,j, i \neq j \\ & q_i \leq u_i \leq Q & \forall i \\ & x_{ij} \in \{0,1\} & \forall i,j \end{array} $$# #random instance
# n = 10
# m = 4
# seed = 1
# random.seed(seed)
# n =10
# x = dict([(i,100*random.random()) for i in range(n)])
# y = dict([(i,100*random.random()) for i in range(n)])
# c, q = { }, {}
# pos ={}
# for i in range(n):
# q[i] = random.randint(10,20)
# pos[i] = x[i],y[i]
# for j in range(n):
# if i !=j:
# c[i,j] = distance(x[i],y[i],x[j],y[j])
folder = "../data/cvrp/"
#file_name = "E-n22-k4.vrp"
#file_name = "E-n30-k3.vrp"
file_name = "E-n51-k5.vrp"
m = 5 #運搬車の台数；ベンチマーク問題例のk?の?を代入
f = open(folder + file_name)
data = f.readlines()
f.close()
n = int(data[3].split()[-1])
Q = int(data[5].split()[-1])
print("n=",n, "Q=",Q)
x, y, q = {}, {}, {}
pos = {}
for i, row in enumerate(data[7 : 7 + n]):
id_no, x[i], y[i] = list(map(int, row.split()))
pos[i] = x[i], y[i]
for i, row in enumerate(data[8 + n : 8 + 2 * n]):
id_no, q[i] = list(map(int, row.split()))
c = {}
for i in range(n):
for j in range(n):
if i != j:
c[i, j] = int(distance(x[i], y[i], x[j], y[j]))
model = Model("avrp - mtz")
x, u = {}, {}
for i in range(n):
u[i] = model.addVar(lb=q[i], ub=Q, vtype="C", name="u(%s)" % i)
for j in range(n):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n):
model.addConstr(quicksum(x[i, j] for j in range(n) if j != i) == 1, "Out(%s)" % i)
model.addConstr(quicksum(x[j, i] for j in range(n) if j != i) == 1, "In(%s)" % i)
# depot
model.addConstr(quicksum(x[0, j] for j in range(1, n) if j != i) == m, "Out(%s)" % 0)
model.addConstr(quicksum(x[j, 0] for j in range(1, n) if j != i) == m, "In(%s)" % 0)
for i in range(n):
for j in range(1, n):
if i != j:
model.addConstr(
u[i] - u[j] + Q * x[i, j] <= Q - q[j], "MTZ(%s,%s)" % (i, j)
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.optimize()
G = nx.Graph()
for (i, j) in x:
if x[i, j].X > 0.001:
# print(i,j)
G.add_edge(i, j)
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
時間枠付き配送計画問題
ポテンシャルを用いた定式化は，時間枠付き配送計画問題に対しても使える．この場合には，ポテンシャルは時刻になる．詳細は，時間枠付き巡回セールスマン問題を参照されたい．
非対称な移動費用をもつ容量制約付きの配送計画問題に対してポテンシャル定式化を拡張した定式化は，以下のように記述できる．
$$ \begin{array}{lll} minimize & \sum_{i,j} c_{ij} x_{ij} & \\ s.t. & \sum_{j} x_{ji} = \sum_{j} x_{ij} = 1 & \forall i \\ & \sum_{j} x_{0j} = \sum_{j} x_{j0} = m & \\ & u_i - u_j + Q x_{ij} \leq Q -q_j & \forall i,j, i \neq j \\ & t_i - t_j + M x_{ij} \leq M -c_{ij} & \forall i,j, i \neq j \\ & q_i \leq u_i \leq Q & \forall i \\ & e_i \leq t_i \leq \ell_i & \forall i \\ & x_{ij} \in \{0,1\} & \forall i,j \end{array} $$時間枠付き配送計画問題に対するベンチマーク問題例は，以下のサイトからダウンロードできる．
https://people.idsia.ch/~luca/macs-vrptw/solutions/welcome.htm
folder = "../data/vrptw/"
#file_name = "RC208.txt"
file_name = "c101.txt"
f = open(folder + file_name)
data = f.readlines()
f.close()
m, Q = list(map(int, data[4].split()))
print("m=",m, "Q=", Q)
x, y, q, e, l, s = {}, {}, {}, {}, {}, {}
pos = {}
for i, row in enumerate(data[9:]):
id_num, x[i], y[i], q[i], e[i], l[i], s[i] = list(map(int, row.split()))
pos[i] = x[i], y[i]
n = len(x)
c = {}
for i in range(n):
for j in range(n):
if i != j:
c[i, j] = distance(x[i], y[i], x[j], y[j]) + s[i] #s[i]は作業時間
model = Model("vrptw - mtz")
x, u, t = {}, {}, {}
for i in range(n):
u[i] = model.addVar(lb=q[i], ub=Q, vtype="C", name="u(%s)" % i)
t[i] = model.addVar(lb=e[i], ub=l[i], vtype="C", name="t(%s)" % i)
for j in range(n):
if i != j:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in range(1, n):
model.addConstr(quicksum(x[i, j] for j in range(n) if j != i) == 1, "Out(%s)" % i)
model.addConstr(quicksum(x[j, i] for j in range(n) if j != i) == 1, "In(%s)" % i)
# depot
model.addConstr(quicksum(x[0, j] for j in range(1, n) if j != i) == m, "Out(%s)" % 0)
model.addConstr(quicksum(x[j, 0] for j in range(1, n) if j != i) == m, "In(%s)" % 0)
for i in range(n):
for j in range(1, n):
if i != j:
M = max(l[i] + c[i, j] - e[j], 0)
model.addConstr(
u[i] - u[j] + Q * x[i, j] <= Q - q[j], "MTZ(%s,%s)" % (i, j)
)
model.addConstr(
t[i] - t[j] + M * x[i, j] <= M - c[i, j], "MTZ2(%s,%s)" % (i, j)
)
model.setObjective(quicksum(c[i, j] * x[i, j] for (i, j) in x), GRB.MINIMIZE)
model.optimize()
G = nx.Graph()
for (i, j) in x:
if x[i, j].X > 0.001:
G.add_edge(i, j)
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
時間枠付き配送計画問題に対するメタヒューリスティクスの設計方法の基本原理
時間枠付き配送計画問題に対する近傍探索ベースのメタヒューリスティクスを設計するためには，補助配列の利用が有効である．
パス $P$ に対して，以下の量を補助配列に保管しておく．
- 総時間 $TT(P)$ (total time)
- 最後の点の最早到着時刻 $ED(P)$ (earliest departure)
- 最初の点の最遅発時刻 $LD(P)$ (latest departure)：その時刻より遅く出ると，パス内のどこかの点で時間枠を逸脱する．
パス $P$ （最終地点を $i$）に点 $j$ （時間枠 $[E_j, L_j]$）を追加する際には， 以下のように更新する． 点 $i,j$ 間の移動時間を $T_{ij}$，点 $i$ の作業時間を $S_i$ とする．
$$ TT := TT + T_{ij} + S_j $$$$ ED := \max \{ ED + T_{ij}, E_j \} + S_j $$$$ LD := \min \{ LD, L_j-TT-T_{ij} \} $$発時刻から $T_{max}$ 時間以内であるための条件は， 最後の地点（デポ）への最早到着時刻が，最遅発時刻に $T_{max}$ を加えた時刻以下であれば良いので，以下の式で確認できる．
$$ ED \leq LD + T_{max} $$トレーラー型配送計画問題（集合被覆アプローチ）
1ルートあたりの配送件数が少ない容量制約付き配送計画問題は，トレーラー型配送計画問題とよばれ， ルートをあらかじめ列挙しておくことによって集合被覆（分割）問題に帰着させることができる．
実行可能なルートの集合を $R$ と書く． ここで実行可能とは，ルートに含まれている顧客を，与えられたすべての条件を満たして巡回可能なことを指す． $a_{ir}$ を顧客 $i$ がルート $r$ に含まれているとき $1$，それ以外のとき $0$ であるパラメータとする． 容量制約付き配送計画問題の場合には，実行可能なルート $r(\in R)$ は以下の式を満たす． $$ \sum_{i \in V_0} a_{ij} q_i \leq Q $$
ルート $r \in R$ に付随する費用を $c_r$ と書く．これは巡回費用と，運搬車の固定費用の和になる． これは，ルートに含まれる顧客の集合に対して最小費用の訪問順を求めることによって決められる． この問題は，デポとルートに含まれる顧客を合わせた点集合に対する巡回セールスマン問題（時間枠が付いている場合には，時間枠付き巡回セールスマン問題）となる． 一般には，この問題は難しいが，ここでは1つのルートに含まれる顧客数が少ないことを想定している （それが集合被覆問題を用いた解法を適用する際の条件となる）ので，比較的容易に解くことができる． 通常は，1つのルートに含まれる顧客数が $5$ 以下なので，全列挙や動的最適化を用いれば良い．
以下の $0$-$1$ 変数を用いる． $$ x_r = \left\{ \begin{array}{ll} 1 & ルート r が最適ルートに含まれる \\ 0 & それ以外 \end{array} \right. $$
定式化は，以下のように書ける．
$$ \begin{array}{ll} minimize & \sum_{r \in R} c_r x_r & \\ s.t. & \sum_{r \in R} a_{ir} x_r \geq 1 & \forall i \in V_0 \\ & x_r \in \{0,1\} & \forall r \in R \end{array} $$上の問題は，制約式の係数が $0$ または $1$ で，右辺定数がすべて $1$で あるという特殊構造をもつ $0$-$1$整数最適化問題であり，一般に集合被覆問題（set covering problem）とよばれる． また，制約がすべて等号の場合を集合分割問題(set partitioning problem)とよぶ，
運搬車の台数を $m$ 台に固定したときの定式化は，以下のように書ける．
$$ \begin{array}{ll} minimize & \sum_{r \in R} c_r x_r & \\ s.t. & \sum_{r \in R} a_{ir} x_r \geq 1 & \forall i \in V_0 \\ & \sum_{r \in R} x_r = m & \\ & x_r \in \{0,1\} & \forall r \in R \end{array} $$この問題は，数理最適化ソルバーで比較的簡単に解くことができる．
以下では，ランダムな問題例に対して，ルートに含まれる点数が[LB,UB]の部分集合を列挙し，集合被覆アプローチを適用する．
random.seed(1)
n = 50
xx = dict([(i, 100 * random.random()) for i in range(n)])
yy = dict([(i, 100 * random.random()) for i in range(n)])
G = nx.Graph()
pos = {}
for i in range(n):
pos[i] = xx[i], yy[i]
G.add_node(i)
V = G.nodes()
c, q = {}, {}
Q = 100
for i in V:
q[i] = random.randint(15, 65)
for j in V:
c[i, j] = distance(xx[i], yy[i], xx[j], yy[j])
G.add_edge(i,j)
q[0] = 0
def generate_routes(V,c,q,Q,LB=1,UB=1,fixed_cost=0):
"""
ルート生成：点数が[LB,UB]の部分集合を列挙して最適順回路を求める
"""
def findsubsets(s, n):
return [set(i) for i in itertools.combinations(s, n)]
n = len(V)
V0 = set(V) - {0}
route = {}
for i in V0:
route[i] = []
cost = {}
visit = {}
r = 0
for n in range(LB, UB + 1):
for S in findsubsets(V0, n):
if sum(q[i] for i in S) > Q:
continue
SS = list(S)
SS.insert(0, 0)
opt_val, tour = tsp.tspdp(len(SS), c, SS)
cost[r] = opt_val + fixed_cost
visit[r] = tour
for i in tour[1:-1]:
route[i].append(r)
r += 1
num_of_routes = r
return cost, route, visit, num_of_routes
cost, route, visit, num_of_routes = generate_routes(V,c,q,Q,LB=1,UB=3,fixed_cost=0)
V0 = set(V) - {0}
model = Model("scp")
x = {}
for r in cost:
x[r] = model.addVar(vtype="B", name=f"x[{r}]")
model.update()
constr = {}
for i in V0:
constr[i] = model.addConstr(quicksum(x[r] for r in route[i]) >= 1, name=f"cover({i})")
model.setObjective(quicksum(cost[r] * x[r] for r in cost), GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
edges = []
for r in x:
if x[r].X > 0.001:
for idx, i in enumerate(visit[r][:-1]):
edges.append((i, visit[r][idx + 1]))
edges.append((visit[r][-1], 0))
G.remove_edges_from(G.edges())
nx.draw(G, pos=pos, with_labels=True, node_color="y")
nx.draw_networkx_edges(G, pos=pos, edgelist=edges, edge_color="blue", width=1);
folder = "../data/cvrp/"
file_name = "E-n22-k4.vrp"
#file_name = "E-n30-k3.vrp"
#file_name = "E-n51-k5.vrp"
m = 4
f = open(folder + file_name)
data = f.readlines()
f.close()
n = int(data[3].split()[-1])
Q = int(data[5].split()[-1])
print("n=",n, "Q=",Q)
x, y, q = {}, {}, {}
pos = {}
for i, row in enumerate(data[7 : 7 + n]):
id_no, x[i], y[i] = list(map(int, row.split()))
pos[i] = x[i], y[i]
for i, row in enumerate(data[8 + n : 8 + 2 * n]):
id_no, q[i] = list(map(int, row.split()))
c = {}
for i in range(n):
for j in range(n):
c[i, j] = int(distance(x[i], y[i], x[j], y[j]))
G = nx.Graph()
for i in range(n):
for j in range(n):
if i!=j:
G.add_edge(i,j)
V = G.nodes()
def beam_search(bests, num_bests, layer):
"""
ビームサーチ（広がり優先）
"""
all_empty_flag = True
bests[layer+1] = []
for (total_prize, travel_time, total_q, a_path) in bests[layer]:
i = a_path[-1]
path = a_path[:]
vals = []
for j in G.neighbors(i):
if j not in path:
reduced_cost = total_prize + c[i,j] -c[i,source] + c[j,source]-prize[j]
if travel_time + c[i,j] + c[j,source]>Tmax or total_q+q[j]>Q: #時間と容量をcheck
pass
else:
vals.append( (reduced_cost,j) )
if len(vals)==0:
#print("empty",i)
pass
else:
all_empty_flag = False
vals.sort(reverse=False)
for reduced_cost,j in vals[:width]:
new_travel_time = travel_time+c[i,j]
bests[layer+1].append( (reduced_cost,new_travel_time, total_q+q[j], path+[j]) )
if all_empty_flag:
return bests
bests[layer+1].sort(reverse=False)
bests[layer+1] = bests[layer+1][:num_bests]
bests = beam_search(bests, num_bests, layer+1)
return bests
# 初期ルート（ピストン輸送）生成
cost, route, visit, num_of_routes = generate_routes(V,c,q,Q,LB=1,UB=1,fixed_cost=0)
model = Model("scp-column generation")
V0 = set(V) - {0}
x = {}
for r in cost:
x[r] = model.addVar(vtype="B", name=f"x[{r}]")
model.update()
constr = {}
for i in V0:
constr[i] = model.addConstr(quicksum(x[r] for r in route[i]) >= 1, name=f"cover({i})")
#model.addConstr(quicksum(x[r] for r in cost) == m, name="num_of_vehicles") #ここで入れると実行不能になるので，列を生成した後に入れる
model.setObjective(quicksum(cost[r] * x[r] for r in cost), GRB.MINIMIZE)
model.Params.OutputFlag = 0
model.update()
while 1:
relax = model.relax()
# relax.optimize(solver = PULP_CBC_CMD(mip=False, presolve=False)) #mypulpの場合ここを生かす
relax.optimize()
dual = {}
for i, con in enumerate(relax.getConstrs()):
dual[i + 1] = con.Pi
dual[0] = 0
source = 0
Tmax = 1000000
prize = dual.copy()
#Beam Search
path = [source]
bests = [ (-prize[source], 0, 0, [source]) ]
width = 200
num_bests = 200
bests = {0: [(-prize[source],0,0,[source])] }
bests_ = beam_search(bests, num_bests, 0)
bests =[]
for k in bests_:
bests.extend( bests_[k] )
bests.sort(reverse=False)
print("min_cost", bests[0])
if bests[0][0]>=-0.000001:
break
for (total_prize, travel_time, slack, tour) in bests:
if total_prize>=-0.000001:
break
cost_ = travel_time+c[tour[-1],source]
col = Column()
for i, cust in enumerate(tour):
if i != 0:
col.addTerms(1, constr[cust])
num_of_routes += 1 # variable index
visit[num_of_routes] = tour
x[num_of_routes] = model.addVar(obj=cost_, vtype="B", column=col)
model.update()
print(num_of_routes)
model.addConstr(quicksum(x[r] for r in x) == m, name="num_of_vehicles")
model.optimize()
print(model.ObjVal)
edges = []
for r in x:
if x[r].X > 0.001:
print(r, x[r].X, visit[r])
for idx, i in enumerate(visit[r][:-1]):
edges.append((i, visit[r][idx + 1]))
edges.append((visit[r][-1], 0))
G.remove_edges_from(G.edges())
nx.draw(G, pos=pos, with_labels=True, node_color="y")
nx.draw_networkx_edges(G, pos=pos, edgelist=edges, edge_color="blue", width=5);
分割配送計画問題
デポからの配送もしくは収集だけを考えた配送計画において，荷の分割を考慮することを考える． ただし，荷を分割することは入庫作業の繁雑さを招くため，固定費用がかかるものとし， さらに，実務においては，荷ごとに分割数の上限が設定される場合もある．
実務においては，そもそも1つの顧客の需要量が，運搬車の積載容量を超えている場合がある． このような場合には，荷の分割は必須になる．最も簡単な方法は，できるだけ積載容量の上限で配送するように分割することである． ただし，これは簡便法であり， 移動費用が三角不等式 ($c_{ik} + c_{kj} \leq c_{ij}$）を満たしていても， 必ずしも最適解になるとは限らない．
特殊な例だが，最適解にならない例を示しておこう． $Q=6, n=7$ で，顧客１の需要が $2Q = 12 $で，他の顧客の需要が$1$であるとする． 移動費用は，顧客1と他の顧客間が $1$，デポ $0$ と1以外の顧客間が$1$であり，他はすべて $2$ とする． このとき，積載量の上限で分割するヒューリスティクスの場合には，顧客１へは２つの往復ルートで配送をし，他の6つの顧客へは1つのルートで巡回するので，総費用は $4 \times 2 +2+ 2\times 5 = 20$ となる． 一方，最適な分割配送は，顧客1へ $4$ ずつ3回に分けて運び，ルートは「デポ・他の顧客・顧客1・他の顧客・デポ」の順とする． この解の費用は $4 \times 3=12$ であり， ヒューリスティクスの近似比率は $20/12$ となる．
簡単な解析から分かるように，分割配送の効果が大きいのは，需要量 $q$ が，運搬車の容量 $Q$ の半分より少し大きい場合である． 移動費用を無視して，容量だけを考えた場合にはビンパッキング問題になるので，ビンパッキング問題の最適値が，分割配送をした場合の運搬車台数 $$ \left\lceil \sum_{i} q_i /Q \right\rceil $$ と比べて，大きい場合には分割配送が有利になる． たとえば， $100$ の顧客が同じ位置にいて，需要量 $q$ がすべて $51$， 容量が $100$ の場合には，ビンパッキング問題の最適値（すべてピストン輸送に相当）は $100$ になるが， 分割配送をすれば2つの顧客を巡回するルートと，すべてのルートを一巡するルートになるので，$51$台の運搬車で処理できる．
荷を分割して配送するのが有利になるのは，荷量が比較的大きい場合であると想定される． このような場合には，運搬車が経由する顧客数は，比較的少ないので， ルートをあらかじめ生成しておく方法（ルート生成法）が有効になる．
集合:
- $R$: デポから出発した運搬車は，幾つかの荷を運び，再びデポに戻るものとする．これをルートとよぶ． ここでは，あらかじめ良いルートの候補を列挙しておくものとし，ルートの集合を $R$ とする．
- $V_0$: 顧客の集合
パラメータ:
- $C_r$: ルート $r (\in R)$ を使用したときにかかる費用（移動費用と固定費用の和）
変数:
- $x_{r}$: ルート $r$ が最適ルートに含まれるとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数．
- $y_{ir}$: 顧客 $i$ の需要をルート $r$ の運搬車によって運んだ量を表す実数変数
上の記号を用いると，定式化は以下のように書ける．
$$ \begin{array}{lll} minimize & \sum_{r \in R} C_r x_r & \\ s.t. & \sum_{r \in R} a_{ir} x_{r} \geq 1 & i \in V_0 \\ & \sum_{r \in R} y_{i r} = q_i & i \in V_0 \\ & y_{ir} \leq q_i x_r & i \in V_0, r \in R \\ & \sum_{i} y_{ir} \leq Q x_r & r \in R \\ & y_{ir} \geq 0 & i, r \in R \\ & x_r \in \{0,1\} & r \in R \end{array} $$上の定式化では，顧客を巡回しても，何も運ばないという解が出る可能性がある． そのような場合には，何も運ばない顧客をルートから顧客を除くことによって，より効率的な解に変換することができる．
random.seed(1)
n = 50
xx = dict([(i, 100 * random.random()) for i in range(n)])
yy = dict([(i, 100 * random.random()) for i in range(n)])
G = nx.Graph()
pos = {}
for i in range(n):
pos[i] = xx[i], yy[i]
G.add_node(i)
V = G.nodes()
c, q = {}, {}
Q = 100
for i in V:
q[i] = random.randint(15, 65)
for j in V:
c[i, j] = distance(xx[i], yy[i], xx[j], yy[j])
G.add_edge(i,j)
q[0] = 0
# Qを超えている部分集合も列挙
route = {}
for i in V0:
route[i] = []
cost = {}
visit = {}
r = 0
LB = 1
UB = 2
for n in range(LB, UB + 1):
for S in findsubsets(V0, n):
SS = list(S)
SS.insert(0, 0)
opt_val, tour = tsp.tspdp(len(SS), c, SS)
cost[r] = opt_val + fixed_cost
visit[r] = tour
for i in tour[1:-1]:
route[i].append(r)
r += 1
# すべての点を通過するルート（ガーベッジコレクション）を追加
cc = {(i, j): int(c[i, j] * 10000) for (i, j) in c}
opt_val, tour, G = tsp.tsplk(len(V), cc)
cost[r] = opt_val / 10000
visit[r] = tour
for i in tour[1:-1]:
route[i].append(r)
r += 1
print("Number of routes=",r)
model = Model("svrp")
x, y = {}, {}
for r in cost:
x[r] = model.addVar(vtype="B", name=f"x[{r}]")
for i in V0:
y[i, r] = model.addVar(vtype="C", name=f"y[{i},{r}]")
model.update()
customer_constr = {}
for i in V0:
customer_constr[i] = model.addConstr(
quicksum(x[r] for r in route[i]) >= 1, name=f"Customer[{i}]"
)
demand_constr = {}
for i in V0:
demand_constr[i] = model.addConstr(
quicksum(y[i, r] for r in route[i]) == q[i], name=f"Demand[{i}]"
)
for i in V0:
for r in route[i]:
model.addConstr(y[i, r] <= q[i] * x[r])
for r in x:
model.addConstr(quicksum(y[i, r] for i in visit[r][1:-1]) <= Q * x[r])
model.setObjective(quicksum(cost[r] * x[r] for r in cost), GRB.MINIMIZE)
model.optimize()
print(model.ObjVal)
何も運ばない顧客をルートから顧客を除く．
changed = set([])
for i in V0:
for r in route[i]:
if x[r].X > 0.001 and y[i, r].X <= 0.001:
print("delete", i, "from", r)
visit[r].remove(i)
changed.add(r)
for r in changed:
opt_val, tour = tspdp(len(visit[r]), c, list(visit[r]))
cost[r] = opt_val
visit[r] = tour
解を描画する．
edges = []
for r in x:
if x[r].X > 0.001:
#print(r, cost[r], x[r].X, visit[r])
for idx, i in enumerate(visit[r][:-1]):
edges.append((i, visit[r][idx + 1]))
edges.append((visit[r][-1], 0))
nx.draw(G, pos=pos, with_labels=True, node_color="y")
nx.draw_networkx_edges(G, pos=pos, edgelist=edges, edge_color="blue", width=1);
分割配送計画問題の別の定式化
実務においては，需要の分割数の上限を設定したいことがある．ここでは，そのような場合の定式化を示す．
集合:
- $L$: デポから顧客への輸送の要求を荷とよぶ．荷の集合を $L$ と定義する
- $R$: デポから出発した運搬車は，幾つかの荷を運び，再びデポに戻るものとする．これをルートとよぶ． ここでは，あらかじめ良いルートの候補を列挙しておくものとし，ルートの集合を $R$ とする
- $L_r$: ルートによって運ぶことができる荷の集合を，ルート内の荷の集合とよぶ． ルート $r(\in R)$ 内の荷の集合を $L_r$ と記す
パラメータ:
- $C_r$: ルート $r (\in R)$ を使用したときにかかる費用
- $F_{\ell}$: 荷 $\ell (\in L)$ の荷受け作業に要する費用． $F_{\ell}$ は分割に対するペナルティを表すので，$F_{\ell}$ が大きいほど，荷を分割しない解を得ることができる
- $U_{\ell}$ : 荷 $\ell (\in L)$ の訪問回数の上限． たとえば，$U_{\ell}=1$ の場合には，高々 $1$ 回の訪問を許すので，荷の分割はできないことを表し， $U_{\ell}=2$ の場合には，荷の分割を $1$ 回だけ許すことを表す
- $D_{\ell}$: 荷 $\ell (\in L)$ の運ぶ量（需要量）
- $M_{r}$: ルート $r (\in R)$ に対応する運搬車の積載重量の上限． ルート（に割り振られた）運搬車の容量を表す
変数:
- $x_{\ell r}$: 荷 $\ell (\in L)$ をルート $r (\in R)$ によって運んだ量を表す実数変数
- $y_{\ell r}$: 荷 $\ell (\in L)$ をルート $r (\in R)$ によって運んだとき $1$， それ以外のとき $0$ を表す $0$-$1$ 変数
- $z_{r}$: ルート $r (\in R)$ を使用したとき $1$，それ以外のとき $0$ を表す $0$-$1$ 変数
上の記号を用いると，定式化は以下のように書ける．
$$ \begin{array}{lll} minimize & \sum_{r \in R} C_r z_r +\sum_{r \in R} \sum_{\ell \in L_r} F_{\ell r} y_{\ell r} & \\ s.t. & \sum_{r \in R: \ell \in L_r} x_{\ell r} =D_{\ell} & \ell \in L \\ & x_{\ell r} \leq \min \{D_{\ell},M_r\} y_{\ell r} & \ell \in L_r, r \in R \\ & \sum_{r \in R: \ell \in L_r} y_{\ell r} \leq U_{\ell} & \ell \in L \\ & \sum_{\ell \in L_r} x_{\ell r} \leq M_r z_r & r \in R \\ & x_{\ell r} \geq 0 & \ell \in L_r, r \in R \\ & y_{\ell r} \in \{0,1\} & \ell \in L_r, r \in R \\ & z_r \in \{0,1\} & r \in R \end{array} $$巡回セールスマン型配送計画問題 （ルート先・クラスター後法）
ルート先・クラスター後法（route-first/cluster-second method）とは，はじめにすべての点（顧客およびデポ）を通過する 巡回路を（たとえば巡回セールスマン問題を解くことによって）作成し，その後でそれをクラスターに分けることによって運搬車のルートを生成する解法の総称である．
積載量制約付きの配送計画問題に対して，すべての点を通過する巡回路を，その巡回路の順番を崩さないように「最適」に分割する方法ができる．
点（顧客およびデポ）の集合 $N$ をちょうど1回ずつ通過する巡回路を表す順列を $\rho$とする． $\rho (i)$ は $i$ 番目に通過する点の番号であり， $\rho(0)$ はデポ $(0)$ である． $C_{ij}$ を $\sum_{k=i+1}^j q_{\rho(k)} \leq Q$ のとき $i+1$ 番目から $j$ 番目の顧客を $\rho$ で定義される順に経由したときのルートの費用とし， それ以外のとき $\infty$ と定義する． すなわち $$ C_{ij} =\left\{ \begin{array}{c l} ルート (0, \rho(i+1), \cdots, \rho(j),0) の費用 & \sum_{k=i+1}^j q_{\rho(k)} \leq Q のとき \\ \infty & それ以外 \end{array} \right. $$ である．
$C_{ij}$ を枝の費用としたとき，点 $0$ から $n=|N|$ までの（有向閉路をもたないグラフ上での）最短路は，動的最適化で計算できる． 最短路に対応する巡回路 $\rho$ が，最適な分割になる．
$j$番目の点までの最適値を $F_j$ とする． $F_0 = 0$ の初期条件の下で，以下の再帰方程式によって最適値を得ることができる．
$$ F_j = \min \{ F_i +C_{ij} \} \ \ j=1,2,\ldots,n $$cost = {
(0, 1): 20,
(0, 2): 25,
(0, 3): 30,
(0, 4): 40,
(0, 5): 35,
(1, 2): 10,
(2, 3): 30,
(3, 4): 25,
(4, 5): 15,
(1, 0): 20,
(2, 0): 25,
(3, 0): 30,
(4, 0): 40,
(5, 0): 35,
}
demand = {1: 5, 2: 4, 3: 4, 4: 2, 5: 7}
Q = 10
n = len(demand)
V = defaultdict(lambda: np.inf)
prev = defaultdict(int)
V[0] = 0
for i in range(1, n + 1):
c = V[i - 1] + cost[0, i]
q = demand[i]
if c + cost[i, 0] < V[i]:
V[i] = c + cost[i, 0]
prev[i] = i - 1
for j in range(i + 1, n + 1):
c += cost[j - 1, j]
q += demand[j]
if q > Q:
break
if c + cost[j, 0] < V[j]:
V[j] = c + cost[j, 0]
prev[j] = i - 1
G = nx.DiGraph()
G.add_nodes_from(list(range(n + 1)))
now = n
while now > 0:
p = prev[now]
G.add_edge(0, p + 1)
for i in range(p + 1, now):
G.add_edge(i, i + 1)
G.add_edge(now, 0)
now = p
pos = nx.circular_layout(G)
nx.draw(G, pos=pos, node_size=1000, with_labels=True, node_color="yellow")
plt.show()
空間充填曲線法
ここで述べる空間充填曲線法 (spacefilling curve method) も ルート先・クラスター後法の一種であるが，手軽に利用できるという利点をもつ．
基本となるのは直線を平面全体に移す曲線（空間充填曲線: spacefilling curve， 怪物曲線: monster curve）である．
- 配送先の $x,y$ 座標を地図から読みとる．
- 各地点(デポを含む)ごとに空間充填曲線上での位置（逆写像） $\theta$ の値を計算する．
- $\theta$ の値を小さい順に並べる．
- デポから順番に配送先を回り運搬車の容量を超えたらデポに戻る．
より簡便な方法としては $x,y,\theta$ を記入した顧客カード（名刺大の厚紙） を $\theta$ の小さい順に箱（日本ではあまり見かけないが，ローロデックスボックスとよばれる名刺入れが最適である） に入れておき，ほぼ均等になるように 運搬車台数分に分割して運転手に渡せば良い．運転手は渡されたカードの順番に配送先を巡回すれば 計算機を使わずに比較的良いルートを見つけることができる．
空間充填曲線上の値を計算する関数sfcを以下に示す．
引数：
- x: x座標
- y: y座標
- maxInput: 座標の最大値（最小値は0と仮定）；既定値は100
返値：
- 空間充填曲線（シェルピンスキー曲線）上の値
random.seed(12)
n = 100
x = {i: 100 * random.random() for i in range(1,n+1)}
y = {i: 100 * random.random() for i in range(1,n+1)}
x[0], y[0] = 50, 50 #depot
theta = []
for i in range(n+1):
theta.append((sfc(x[i], y[i]), i))
theta.sort()
tour = []
for _, i in theta:
tour.append(i)
G = nx.Graph()
for idx, i in enumerate(tour[:-1]):
G.add_edge(i, tour[idx + 1])
G.add_edge(tour[-1], tour[0])
pos = {i: (x[i], y[i]) for i in range(n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
上の空間充填曲線法で得た順回路に対して，ルート先・クラスター後法を適用する．
demand = {i:random.randint(1,10) for i in range(1,n+1) }
cost = { (i,j): distance(x[i], y[i], x[j], y[j]) for i in range(n+1) for j in range(n+1) }
Q = 50
V = defaultdict(lambda: np.inf)
prev = defaultdict(int)
V[0] = 0
#0(depot)を最初にした巡回順tour0を求める
idx = tour.index(0)
tour0 = tour[idx:] + tour[:idx]
for i in range(1,n+1):
cust_i = tour0[i]
c = V[i - 1] + cost[0, cust_i]
q = demand[cust_i]
if c + cost[cust_i, 0] < V[i]:
V[i] = c + cost[cust_i, 0]
prev[i] = i - 1
for j in range(i + 1, n + 1):
cust_j = tour0[j]
c += cost[tour0[j - 1], cust_j]
q += demand[cust_j]
if q > Q:
break
if c + cost[cust_j, 0] < V[j]:
V[j] = c + cost[cust_j, 0]
prev[j] = i - 1
G = nx.Graph()
G.add_nodes_from(list(range(n + 1)))
now = n
while now > 0:
p = prev[now]
G.add_edge(0, tour0[p + 1])
for i in range(p + 1, now):
G.add_edge(tour0[i], tour0[i + 1])
G.add_edge(tour0[now], 0)
now = p
pos = {i: (x[i], y[i]) for i in range(n+1)}
nx.draw(G, pos=pos, node_size=1000 / n + 10, with_labels=False, node_color="blue")
plt.show()
運搬車スケジューリング問題
配送計画問題が運搬車経路問題 (vehicle routing problem)と呼ばれるのに対して，ここでは運搬車のスケジューリングを考える． 「スケジューリング」は，時間に焦点を当てたモデリングを行うことを意味する．特に，顧客の訪問時刻が，決められた時刻である場合は，時刻が決められていない（もしくは時間枠として与えられている）配送計画問題と比べて，自由度が減るので，求解が容易になる．
運搬車の台数を最小にする問題は，乗務員スケジューリングの制約なし版とも考えられる．
- $n$: タスクの数
- $s_i, f_i$ : タスク $i$ の開始時刻と終了時刻
- $x_{ij}$: タスク $i$ の直後にタスク $j$ を処理しとき $1$、それ以外のとき $0$
タスクを点、タスク間の移動を枝としたネットワーク $D$ を作る．また，ダミーの始点 $0$ と終点 $n+1$ を追加しておく． このグラフは閉路をもたない有向グラフである．
$$ \begin{array}{lll} minimize & \sum_{i,j} x_{0j} & \\ s.t. & \sum_{j} x_{ji} = \sum_{j} x_{ij} & \forall i=1,2,\ldots,n \\ & \sum_{j} x_{ij} = 1 & \forall i=1,2,\ldots,n \end{array} $$このように発着時刻の指定がついている場合には，どのタスクの次にどのタスクの処理が可能であるかを表すグラフを作成する． 以下の例では，$A,B,C,D,E$ の $5$ つのタスクとデポ（とそのコピー）から成る点集合を生成する． タスク $i$ の着時刻に，タスク $i$ から $j$ への移動時間を加えたものが， タスク $j$ の発時刻以下であるとき，タスク $i$ の後にタスク $j$ が処理可能であるとよび，タスク $i$ を表す点から， タスク $j$ を表す点へ枝をはる．また，デポを表す始点（例題では点 $6$ と仮定する）から，各タスクを表す点に枝をはり， 同様に，各タスクを表す点から終点を表す点のコピーへ枝をはる．
D = nx.DiGraph()
D.add_edges_from(
[("A", "C"), ("A", "D"), ("A", "E"), ("B", "C"),
("B", "D"), ("B", "E"), ("C", "E")]
)
このグラフ上で，始点から終点への独立パス（同じ点を通過しないパス）で，すべての点を被覆する（すべてのタスクを表す点をちょうど $1$回づつ 通過する）ものを求める．実は，最小のパスの本数（運搬車の台数）は， Dilworthの定理とよばれるグラフ理論の結果から，同時に処理できないタスクの最大値と一致することが知られている．
この点を被覆する最小本数の独立パスを求める問題は，マッチング問題を用いて解くことができる． まず，タスクを表す点 $i$ を無向枝 $(i^L, i^R)$ に変換する． さらに，もとのグラフで $i,j$ の枝があるとき $(i^R, j^L)$ の枝をはる． このグラフ上で最大（位数）マッチングを求める．マッチングに含まれる枝が，元のグラフでのパスに含まれる枝に対応する． 各タスクを1台の運搬車で処理する場合と比べて，マッチングの位数 $m$ 分だけ台数が減るので，最小のパス（運搬車）の数は $n-m$ になる．
この例では，同時に処理できないのタスクの最大値は $2$ なので，$2$本の独立パスが存在する． これは，$2$台の運搬車で巡回するのが最適であることを示している．
タスク間に移動費用が定義されている場合も同様であり，マッチングのかわりに重み付きマッチングを求めれば良い．
なお．複数デポだと $NP$-困難になる．
G = nx.Graph()
left = []
pos = {}
for idx, i in enumerate(D.nodes()):
G.add_node(str(i) + "L")
pos[str(i) + "L"] = 0, idx
left.append(str(i) + "L")
G.add_node(str(i) + "R")
pos[str(i) + "R"] = 1, idx
for (i, j) in D.edges():
G.add_edge(str(i) + "R", str(j) + "L")
matching = nx.algorithms.matching.max_weight_matching(G, maxcardinality=True)
plt.figure()
nx.draw(G, pos=pos, node_size=300, with_labels=True, node_color="yellow")
nx.draw_networkx_edges(G, pos=pos, edgelist=list(matching), edge_color="red", width=2)
plt.show()
edges = []
for m in matching:
(i, j) = m
if i[-1] == "R":
edges.append((i[:-1], j[:-1]))
else:
edges.append((j[:-1], i[:-1]))
plt.figure()
pos = nx.circular_layout(D)
nx.draw(D, pos=pos, node_size=100, with_labels=True, node_color="yellow")
nx.draw_networkx_edges(D, pos=pos, edgelist=edges, edge_color="red", width=2)
plt.show()
会議室割当問題への応用
時間帯があたえられた会議を，最小の部屋数で割り振る問題も同じ構造をもっている．これは，https://amplify.fixstars.com/demo で8つの部屋に割り振る例題として使われていたが， 最適値は，以下に示すように5である．
schedules = {
"meeting1": ["10:00", "13:00"],
"meeting2": ["10:00", "12:00"],
"meeting3": ["10:00", "11:00"],
"meeting4": ["11:00", "13:00"],
"meeting5": ["11:00", "12:00"],
"meeting6": ["11:00", "15:00"],
"meeting7": ["12:00", "16:00"],
"meeting8": ["12:00", "15:00"],
"meeting9": ["13:00", "15:00"],
"meeting10": ["13:00", "14:00"],
"meeting11": ["14:00", "17:00"],
"meeting12": ["15:00", "19:00"],
"meeting13": ["15:00", "17:00"],
"meeting14": ["15:00", "16:00"],
"meeting15": ["16:00", "18:00"],
"meeting16": ["16:00", "18:00"],
"meeting17": ["17:00", "19:00"],
"meeting18": ["17:00", "18:00"],
"meeting19": ["18:00", "19:00"],
}
# 会議の数
Nm = len(schedules)
# 会議室の数
Nr = 8
# 時刻を時間単位の数値に変換
def time2num(time: str):
h, m = map(float, time.split(":"))
return h + m / 60
st, fi = {}, {}
for i, m in enumerate(schedules):
st[i], fi[i] = time2num(schedules[m][0]), time2num(schedules[m][1])
D = nx.DiGraph()
for i in range(Nm):
D.add_node(str(i))
for i in st:
for j in st:
if fi[i] <= st[j]:
D.add_edge(str(i), str(j))
G = nx.Graph()
for i in D.nodes():
G.add_node(str(i) + "L")
G.add_node(str(i) + "R")
for (i, j) in D.edges():
G.add_edge(str(i) + "R", str(j) + "L")
matching = nx.algorithms.matching.max_weight_matching(G, maxcardinality=True)
edges = []
for m in matching:
(i, j) = m
if i[-1] == "R":
edges.append((i[:-1], j[:-1]))
else:
edges.append((j[:-1], i[:-1]))
plt.figure()
pos = nx.circular_layout(D)
nx.draw(D, pos=pos, node_size=100, with_labels=True, node_color="yellow", width=0.3)
nx.draw_networkx_edges(D, pos=pos, edgelist=edges, edge_color="red", width=3)
plt.show()
D = nx.DiGraph()
D.add_edges_from(edges)
first_nodes = []
for i in D.in_degree():
# print(i)
if i[1] == 0:
first_nodes.append(i[0])
paths = []
for i in first_nodes:
path = [i]
while True:
succ = list(D.successors(i))
if len(succ) == 0:
break
j = succ[0]
path.append(j)
i = j
paths.append(path)
room_assignment = []
for i, path in enumerate(paths):
for j in path:
room_assignment.append((int(j), i))
num_rooms = len(paths)
mtg_names = list(schedules.keys())
room_names = ["Room " + chr(65 + i) for i in range(num_rooms)]
cmap = plt.get_cmap("coolwarm", num_rooms)
colors = [cmap(i) for i in range(num_rooms)]
fig, ax1 = plt.subplots(nrows=1, ncols=1, figsize=(15, 10))
for mtg_idx, room in room_assignment:
mtg_name = mtg_names[mtg_idx]
start = time2num(schedules[mtg_name][0])
end = time2num(schedules[mtg_name][1])
plt.fill_between(
[room + 0.55, room + 1.45],
[start, start],
[end, end],
edgecolor="black",
linewidth=3.0,
facecolor=colors[room],
)
plt.text(
room + 0.6, start + 0.1, f"{schedules[mtg_name][0]}", va="top", fontsize=10
)
plt.text(
room + 1.0, (start + end) * 0.5, mtg_name, ha="center", va="center", fontsize=15
)
ax1.yaxis.grid()
ax1.set_xlim(0.5, len(room_names) + 0.5)
ax1.set_ylim(19.1, 7.9)
ax1.set_xticks(range(1, len(room_names) + 1))
ax1.set_xticklabels(room_names)
ax1.set_ylabel("Time")
ax2 = ax1.twiny().twinx()
ax2.set_xlim(ax1.get_xlim())
ax2.set_ylim(ax1.get_ylim())
ax2.set_xticks(ax1.get_xticks())
ax2.set_xticklabels(room_names)
ax2.set_ylabel("Time")
plt.show()
積み込み・積み降ろし型配送計画問題
古典的な配送計画問題においては，単一のデポ（配送拠点）から顧客への配送，もしくは顧客からデポへの収集のいずれかを仮定していた． 自然な拡張として，デポから顧客，顧客からデポへの荷の移動を同時に考慮することがあげられる．
さらなく拡張として，デポ以外の地点で積み込み，デポ以外の別の地点で積み降ろしをする荷の輸送も考えられる． これを考慮した問題を積み込み・積み降ろし型配送計画問題(pickup-delivery vehicle routing problem)とよぶ． 積み込み地点から積み降ろし地点までの輸送が直送のみの場合には，空輸送最小化問題になるが， 直送でない場合には問題は複雑になる．
実際問題においては，荷の積み合わせや運搬車のタイプ（横開きのウィング型か否か）を考慮したり，乗り合いタクシー問題(dial-a-ride routing problem)への応用の場合には， 積み込みから積み降ろしまでの経過時間なども考慮する必要が出てくる．
複数デポ配送計画問題
実際問題においては，複数のデポを考慮する必要がある場合が多い．より一般的には，各運搬車の出発地点と到着地点を与えて，デポの概念を取り除いたモデルも可能である．
その他にも，法定の休憩時間の確保， 入庫不能な運搬車と顧客の組合せ，複数次元の需要と容量など，実務では様々な条件が付加される．
すべての条件を考慮することは難しいが，以下の諸制約を考慮した配送最適化システムとして METRO (https://www.logopt.com/metro/) が開発されている．
- 複数時間枠制約付き
- 多次元容量非等質運搬車
- 配達・集荷
- 積み込み・積み降ろし
- 複数休憩条件付き
- スキル条件付き
- 優先度付き
- パス型許容
- 複数デポ（運搬車ごとの出発・到着地点）