数学やプログラミングの備忘録

数理最適化, Python, C++をメインに紹介するブログ。

MENU
トップ > 数理最適化 > 配送計画問題をpythonで最適化する

配送計画問題をpythonで最適化する

数理最適化 数理最適化-配送計画

『容量制約付き配送計画問題(Capacitated Vehicle Routeing Problem, CVRP)を定式化し、pythonのpulpで最適化する方法』を紹介します。

目次

容量制約付き配送計画問題(CVRP)の定式化

まずは、容量制約付き配送計画問題(CVRP)の問題設定を紹介します。

N 人の顧客に対して、デポから複数のトラックで荷物を配達したい。 このとき以下の条件を満たし、配送ルートの距離の合計を最小化せよ:
- それぞれの顧客の需要量だけ配達する
- それぞれの顧客には一度だけ訪問できる
- トラックは、デポから出発しデポに帰って来る
- トラックに積める荷物数に制限がある


続いて、容量制約付き配送計画問題(CVRP)をを整数計画問題に定式化してみます。いろいろ文献やWebページを読み漁ると、定式化の方法は様々あります。どんな定式化にしても、部分巡回路除去制約というものが必要になってくるようです。例えば、下記のWebページを参考にさせていただくと、

運搬経路問題（配送最適化問題，Vehicle Routing Problem) をPuLPで解く - Qiita

このページではDFJ(Dantzig-Fulkerson-Johnson)制約により、部分巡回路路を除去しています。同Webページで名前だけ紹介されていますが、他にも部分巡回路除去制約としてMTZ制約というものがあります。

MTZ制約を調べてみて分かったことは、DFJ制約は {2}^{N} 個のの制約で、これって N が大きいと爆発的に数が増えるじゃん！ってことで、その後提案されたのがMTZ制約で、こちらは {n}^{2} 個の制約となります。

という訳で、本記事は、部分巡回路制約としてMTZ制約を採用します。

まずは、与えられる定数等の記号の定義から、

  • N: 顧客数
  • V = {0,1,\ldots,N}: ノード集合(ただし、0はデポを表す)
  • c_{ij}\ (^\forall i,j\in V with i\neq j): ノード i と ノード j の距離
  • C: トラックに積める荷物数
  • d_i\ (i\in V\setminus{0}): 顧客 i の需要量

次に、整数計画問題の変数を定義します:

\begin{align} x_{ij} &= \begin{cases} 1 \text{ if いずれかのトラックがノード $i$ からノード $j$ へ移動する}\\ 0 \text{ そうでないとき } \end{cases}\\ u_i &= (\text{顧客 $i$ へ荷物を配達したトラックがその時点で配達した荷物量}) \end{align}


MTZ制約(Miller-Tucker-Zemlin subtour elimination constraints) を用いて、CVRPを整数計画問題に定式化していきます。

配送計画問題の定式化 \begin{array}{ll} \min_{x} & \displaystyle\sum_{i\in V}\sum_{j\in V, i\neq j} c_{ij} x_{ij},\\ \text{s.t.} & \displaystyle\sum_{i \in V \setminus \{j\}} x_{ij} = 1 \ (^\forall j \in V \setminus \{0\}),\\ & \displaystyle\sum_{j \in V \setminus \{i\}} x_{ij} = 1 \ (^\forall i \in V \setminus \{0\}),\\ &u_i - u_j + C x_{ij} \le C - d_j\\ &\ (^\forall i,j \in V \setminus \{0\}, i\neq j, \text{such that } d_i + d_j \le C),\\ &d_i \le u_i \le C\ (^\forall i\in V \setminus \{0\}),\\ &x_{ij}\in \{0,1\}\ (^\forall i,j\in V, i\neq j),\\ &u_i \in \mathbb{Z}\ (^\forall i\in V \setminus \{0\}). \end{array}


2種類の制約条件 \sum x = 1 がありますが、これらは、必ず顧客には一度のみ立ち寄ることを意味します。その下の制約条件は、MTZ制約と呼ばれるもので、部分巡回路を除去してくれます。 d_i \le u_i \le C はトラックの容量制約となります。

補足ですが、使用できるトラック数に上限がある場合(上限値 k)は、次の制約条件を入れることになります。

\sum_{j\in V} x_{0j} \le k, \sum_{i\in V} x_{i0} \le k


これらは、それぞれ、デポから出発とデポへ到着できるトラック数が k 以下という意味です。

容量制約付き配送計画(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")

出力結果

f:id:cresselia2012:20200830172406p:plain

いい感じにトラック別で経路に色をつけることができました。

ちなみに、インスタンスがあるWebサイトに最適解があるので、別途可視化してみました。

f:id:cresselia2012:20200830172342p:plain

最適解は、うーん、何というか、綺麗ですね。

まとめ

本記事では、『容量制約付き配送計画問題(Capacitated Vehicle Routeing Problem, CVRP)を定式化し、pythonのpulpで最適化する方法』 を紹介しました。pulpは手軽に使えるので、コンパクトなコードで解いてみることができます。

今回は、顧客数31のインスタンスを解いてみました。最適解まで到達はできませんでしたが、以下のような場合なら最適解へ到達するかもしれません:

  • もっと長い時間をかけて解く
  • 商用ソフトウェアを使う

また、配送計画問題はNP困難と呼ばれるクラスに位置づけられ、難しい問題です。なので、もっと規模の小さい問題であれば、厳密解法であるCBC(pulpに同梱されているソルバー)でも解けるかもしれません。

本記事で紹介したpythonコードはGitHubにもまとめておりますので、試してみたい方はどうぞお使いください。

github.com