『容量制約付き配送計画問題(Capacitated Vehicle Routeing Problem, CVRP)を定式化し、pythonのpulpで最適化する方法』を紹介します。
目次
容量制約付き配送計画問題(CVRP)の定式化
まずは、容量制約付き配送計画問題(CVRP)の問題設定を紹介します。
- それぞれの顧客の需要量だけ配達する
- それぞれの顧客には一度だけ訪問できる
- トラックは、デポから出発しデポに帰って来る
- トラックに積める荷物数に制限がある
続いて、容量制約付き配送計画問題(CVRP)をを整数計画問題に定式化してみます。いろいろ文献やWebページを読み漁ると、定式化の方法は様々あります。どんな定式化にしても、部分巡回路除去制約というものが必要になってくるようです。例えば、下記のWebページを参考にさせていただくと、
運搬経路問題(配送最適化問題,Vehicle Routing Problem) をPuLPで解く - Qiita
このページではDFJ(Dantzig-Fulkerson-Johnson)制約により、部分巡回路路を除去しています。同Webページで名前だけ紹介されていますが、他にも部分巡回路除去制約としてMTZ制約というものがあります。
MTZ制約を調べてみて分かったことは、DFJ制約は 個のの制約で、これって
が大きいと爆発的に数が増えるじゃん!ってことで、その後提案されたのがMTZ制約で、こちらは
個の制約となります。
という訳で、本記事は、部分巡回路制約としてMTZ制約を採用します。
まずは、与えられる定数等の記号の定義から、
: 顧客数
: ノード集合(ただし、0はデポを表す)
with
: ノード
と ノード
の距離
: トラックに積める荷物数
: 顧客
の需要量
次に、整数計画問題の変数を定義します:
MTZ制約(Miller-Tucker-Zemlin subtour elimination constraints) を用いて、CVRPを整数計画問題に定式化していきます。
2種類の制約条件 がありますが、これらは、必ず顧客には一度のみ立ち寄ることを意味します。その下の制約条件は、MTZ制約と呼ばれるもので、部分巡回路を除去してくれます。
はトラックの容量制約となります。
補足ですが、使用できるトラック数に上限がある場合(上限値 )は、次の制約条件を入れることになります。
これらは、それぞれ、デポから出発とデポへ到着できるトラック数が 以下という意味です。
容量制約付き配送計画(CVRP)をpythonのpulpで最適化
pulpは、COIN-ORが開発したpythonのモジュールで、最適化問題を定式化することができます。本記事では、pulpに同梱されているCBCというソルバーで解きます。pulpはターミナルでインストールできます。
pip install pulp
まずは、以下のWebサイトで公開されている Augerat et al. Set A というCVRPのインスタンス(顧客数31)を解いてみます。
Vehicle Routing Problem | NEO Research Group
CVRPをpulpで定式化する前に、必要な係数などを定義します:
import math def makeCVRP(): num_nodes = 32 capacity = 100 demand = [ 0, 19, 21, 6, 19,\ 7, 12, 16, 6, 16,\ 8, 14, 21, 16, 3,\ 22, 18, 19, 1, 24,\ 8, 12, 4, 8, 24,\ 24, 2, 20, 15, 2,\ 14, 9 ] coordinate = [] coordinate.append( ( 82, 76 ) ) coordinate.append( ( 96, 44 ) ) coordinate.append( ( 50, 5 ) ) coordinate.append( ( 49, 8 ) ) coordinate.append( ( 13, 7 ) ) coordinate.append( ( 29, 89 ) ) coordinate.append( ( 58, 30 ) ) coordinate.append( ( 84, 39 ) ) coordinate.append( ( 14, 24 ) ) coordinate.append( ( 2, 39 ) ) coordinate.append( ( 3, 82 ) ) coordinate.append( ( 5, 10 ) ) coordinate.append( ( 98, 52 ) ) coordinate.append( ( 84, 25 ) ) coordinate.append( ( 61, 59 ) ) coordinate.append( ( 1, 65 ) ) coordinate.append( ( 88, 51 ) ) coordinate.append( ( 91, 2 ) ) coordinate.append( ( 19, 32 ) ) coordinate.append( ( 93, 3 ) ) coordinate.append( ( 50, 93 ) ) coordinate.append( ( 98, 14 ) ) coordinate.append( ( 5, 42 ) ) coordinate.append( ( 42, 9 ) ) coordinate.append( ( 61, 62 ) ) coordinate.append( ( 9, 97 ) ) coordinate.append( ( 80, 55 ) ) coordinate.append( ( 57, 69 ) ) coordinate.append( ( 23, 15 ) ) coordinate.append( ( 20, 70 ) ) coordinate.append( ( 85, 60 ) ) coordinate.append( ( 98, 5 ) ) def computeDistance( c1, c2 ): return math.sqrt( pow( c2[0] - c1[0], 2 ) + pow( c2[1] - c1[1], 2 ) ) distance = [ [ round(computeDistance( c1, c2 )) for c1 in coordinate ] \ for c2 in coordinate ] return num_nodes, capacity, demand, distance, coordinate # make problem num_nodes, capacity, demand, distance, coordinate = makeCVRP()
変数表
変数名 | 意味 |
---|---|
num_nodes |
ノード数(顧客数+1) |
capacity |
トラックの容量上限 |
demand |
顧客の荷物需要量ベクトル(デポは0とする) |
distance |
ノード間の距離行列 |
coordinate |
ノードの座標(最適化問題を解く際は必要無い,可視化で必要) |
次に、CVRPをpulpで定式化していきます:
import pulp # 最適化問題を定義 problem = pulp.LpProblem( "CVRP", pulp.LpMinimize ) ## 変数を定義 x = [ [ pulp.LpVariable( 'x_{}_{}'.format( i, j ), cat="Binary" ) \ if i != j else None for j in range(num_nodes) ] \ for i in range(num_nodes) ] u = [ pulp.LpVariable( 'u_{}'.format( i ), demand[i], capacity, cat="Integer" ) \ for i in range(1,num_nodes) ] ## 最小化したい関数を定義 problem += pulp.lpSum( distance[i][j] * x[i][j] for i in range(num_nodes) \ for j in range(num_nodes) if i != j ) ## 制約条件を定義 ### \sum x_ij = 1 for j in range(1,num_nodes): problem += pulp.lpSum( x[i][j] for i in range(num_nodes) if i != j ) == 1 for i in range(1,num_nodes): problem += pulp.lpSum( x[i][j] for j in range(num_nodes) if i != j ) == 1 ### MTZ制約 for i in range(1,num_nodes): for j in range(1,num_nodes): if i != j and demand[i] + demand[j] <= capacity: problem += u[i-1] - u[j-1] + capacity * x[i][j] <= capacity - demand[j] # solve threads = 8 maxSeconds = 60 result = problem.solve(pulp.PULP_CBC_CMD(threads=threads, maxSeconds=maxSeconds)) print("objective value = {}".format(pulp.value(problem.objective)))
かなり短いコードで書くことができました。solve()
関数により、CBCで求解スタートとなります。注意点としては、デフォルトはスレッド数1なので、上記のコードのように threads
で指定しスピードアップしています。また、時間制限は maxSeconds
で秒数指定することができますが、実時間で無く、CPU時間の時間制限 なので注意が必要です。
以下の記事でpulpの実行オプションを解説しているので、興味があればご覧ください:
pulpを使いこなすための備忘録 - 数学やプログラミングの備忘録
検証結果と配送ルートの可視化
上記のコードで、CPU時間3万秒を時間制限とし(maxSeconds = 30000
) 解いてみました。
CBCのログを一部記載します:
Result - Stopped on time limit Objective value: 941.00000000 Lower bound: 453.049
意味を調べると、
- 時間制限で終了
- 目的関数値(配達経路の距離の合計)が941の解を発見
- 最適値の下界は453.049
つまり、時間切れで最適解の保証までは完了していません。 実際、このインスタンスを公開していたWebサイトによると、最適値は784です**。
さて、折角なので、結果を可視化してみます。
networkx と matplotlib が必要ですので、インストールします:
pip install networkx matplotlib
先程のpythonコードの続きです:
# pulp(CBC)の結果から辺を定義 edges = [ ( i, j ) for i in range(num_nodes) for j in range(num_nodes) if i != j and pulp.value(x[i][j]) == 1 ] # edgesから経路を取得 paths = [] for i,j in edges: if i == 0: path = [ i, j ] while path[-1] != 0: for v, u in edges: if v == path[-1]: path.append(u) break paths.append(path) import networkx as nx import matplotlib.pyplot as plt # 有向グラフの作成 G = nx.DiGraph() G.add_edges_from( edges ) color = [ "r", "g", "y", "m", "c" ] edge_color = [] # 経路毎に色をセット for i,j in G.edges: for t,path in enumerate(paths): if i in path and j in path: edge_color.append( color[t] ) break assert len(edges) == len(edge_color) # グラフの描画 pos = { i : coordinate[i] for i in range(num_nodes) } fig = plt.figure() nx.draw_networkx( G, pos, edge_color=edge_color, alpha=0.5) # 画像を保存 plt.axis("off") fig.savefig("test.png")
出力結果
いい感じにトラック別で経路に色をつけることができました。
ちなみに、インスタンスがあるWebサイトに最適解があるので、別途可視化してみました。
最適解は、うーん、何というか、綺麗ですね。
まとめ
本記事では、『容量制約付き配送計画問題(Capacitated Vehicle Routeing Problem, CVRP)を定式化し、pythonのpulpで最適化する方法』 を紹介しました。pulpは手軽に使えるので、コンパクトなコードで解いてみることができます。
今回は、顧客数31のインスタンスを解いてみました。最適解まで到達はできませんでしたが、以下のような場合なら最適解へ到達するかもしれません:
- もっと長い時間をかけて解く
- 商用ソフトウェアを使う
また、配送計画問題はNP困難と呼ばれるクラスに位置づけられ、難しい問題です。なので、もっと規模の小さい問題であれば、厳密解法であるCBC(pulpに同梱されているソルバー)でも解けるかもしれません。