第3回:はじめての配送計画の列生成法【ブレインパッドの数理最適化ブログ】

ブレインパッドの社員が「数理最適化技術」に関して連載するこの企画。
第3回は、当社のデータサイエンティストが、与えられた条件を満たしつつトラックの移動距離が最小になるような配送計画の問題を題材に、列生成法をご紹介します。

f:id:brainpad-inc:20201023143041j:plain

こんにちは。アナリティクスサービス部の岡﨑です。 この記事では 【連載】ブレインパッドの数理最適化ブログ(目次) - Platinum Data Blog by BrainPadの第3弾として、配送計画問題を題材にして列生成法をご紹介いたします。数理最適化に若干の心得があるが列生成法を知らない読者*1が「自分でも列生成を試してみたい!」と感じていただくことを目標に、本記事を執筆しています。

はじめに 

ビジネス課題を解くためのアプローチの一つに、混合整数計画問題(以下MIP)があります。MIPの定式化により、意外なほど多くのビジネス課題が表現できます。また定式化のコツも数多く知られています*2。今回は、列生成法というMIPを解くための手法をご紹介します。MIPを複数の問題に分割し定式化する枠組みの一つです。問題の分割により、素朴には解けない大規模な問題が、現実的な時間で現実的な解を得られる場合があるようです。

容量制約付き配送計画問題 

ここからは、具体的な問題を通して、列生成法の説明します。途中までは最適化に関する予備知識がなくとも、理解できるように書いてみたつもりです。記事は長くなりますが、良ければお付き合いください。

問題設定

1つの倉庫から、複数のトラックを用いて、複数の店舗に荷物を配送する問題を考えます。前提となる条件は下記を想定します。

  • トラックは、倉庫から出発し、倉庫へ帰ることとする。
  • 配送する荷物の量は、店舗ごとに決まっており、各店舗の荷物を1台のトラックで運ぶこととする。
  • トラックに積める最大積載容量が定まっている。

本記事では、与えられた条件を満たしながら、トラックの移動距離を最小にするような配送計画をもとめる方法を考えます。

方針

この問題を解く方法は様々あります*3。今回は列生成法を用いた定式化をします。大まかな方針は、「トラック1台の配送ルートを列挙し、その組み合わせで解(=配送計画)を得る」です。

f:id:bp-writer:20201021095227p:plain
配送ルートの組み合わせで、配送計画を作るイメージ

トラックの積載容量に関する制約については説明しませんでしたが、列挙する配送ルートは、訪問先の荷物を運ぶことが可能なルートのみを考えることにします。

定式化

まずは、配送ルート集合が所与の前提で、その組み合わせを決定する問題を考えます。 知りたいことは、どの配送ルートを採用するかなので、それを表す変数x_rを準備します。


\begin{aligned}
x_r = \begin{cases}1 \quad 配送ルートrを使う\\0 \quad 配送ルートrを使わない\end{cases}
\end{aligned}

この変数を用いて、今回の問題を定式化すると、下記の通りになります。


\begin{aligned}
\mathrm{minimize}&\quad \sum_{r\in R}d_r x_r\\
\mathrm{subject\ to}&\quad \sum_{r\in R}a_{ri} x_r \geq 1,\quad i\in \mathrm{S}\setminus\{0\},\\
&\quad x_{r}\in\{0, 1\},\quad r\in R.\\
\end{aligned}

ただし、

  • S = \{0,1,2,\cdots N\}:倉庫及び店舗の集合(0を倉庫とする)
  • R:配送ルートの集合
  • d_r:配送ルートrでのトラックの総移動距離
  • a_{ri}:配送ルートrに店舗iが含まれている場合に1。それ以外は0

定式化内の各式の表していることは下記の通りです。

  • 目的関数:全トラックの総移動距離(=使う配送ルートに関する移動距離の和)を最小化する
  • 制約条件:配送ルート選択は、全ての配送店舗をカバーできるようにする

非常にシンプルな定式化です。

全ルート列挙

ここからは、定式化に現れる配送ルート集合を作成する方法を考えます。配送ルートを列挙すると全部でどのくらいあるでしょうか?粗い見積もりですが、店舗数をNとすると、\sum_{i=1}^N P(N,i)と書けます。かなり多そうです。

店舗数 配送ルート数
3 15
5 325
10 9.86\times 10^{6}
50 8.26\times 10^{64}

店舗数が10でも配送ルートの総数は、約986万通りあります。配送ルートを全列挙するのは困難です。ただ、実際には必ずしも全ての配送ルートを列挙が必要ないと考えることもできます。例えば5店舗の配送計画であれば、配送ルートの総数325通りのルートを考慮せずとも、パソコンなしで人の手で作れそうです。つまり、全配送ルートの大半は、使われる見込みの少ない不自然なルート(=例えば移動距離が大きすぎるようなルート)と予想できます。例えば、下図のルートは配送には使えなさそうです。

f:id:bp-writer:20201020122605p:plain
使われる見込みの少ない配送ルートのイメージ
まとめると、「全配送ルートは列挙するには多すぎる。ほとんどのルートは最終的な解で使われる見込みが低いであろう。」となります。列生成法は、使われる見込みの高い配送ルートを効率的に列挙する枠組みです。

列生成法による配送ルート列挙

ここから、列生成法について説明いたします。*4。以降の説明は、線形計画問題の双対問題に関する予備知識を仮定します。 先ほどの最適化の定式化を再掲します。


\begin{aligned}
\mathrm{minimize}&\quad \sum_{r\in R}d_r x_r\\
\mathrm{subject\ to}&\quad \sum_{r\in R}a_{ri} x_r \geq 1,\quad i\in \mathrm{S}\setminus\{0\},\\
&\quad x_{r}\in\{0, 1\},\quad r\in R.\\
\end{aligned}

配送ルート集合Rが所与で、かつ、定式化の解を求めることが出来れば、解きたかった問題の解が得られます。これまでの議論で、「配送ルートを全て列挙することは困難であるが、不要そうな配送ルートを除外した配送ルートの候補の部分集合から求解しても、元の最適化問題の解が変わらない」ことを予想しました。 ここからは、R_k\subset Rを満たす、良い配送ルート集合R_kを求めることを考えます。

列生成法の方針

まず、R_0なる初期配送ルート集合を考えます。初期集合に新たな配送ルートを追加していき、良い配送ルート集合を作ることを目指します。良い配送ルートを作る素朴な方法は、「新しい配送ルートを作成し、既存の配送ルート集合に追加して、元の最適化問題を解く、という手順を配送計画の解が改善しなくなるまで繰り返す」です。しかしながら、この素朴な方法では、バイナリ変数を含む最適化問題を繰り返し解く必要があり、計算効率が良くありません。

線形緩和問題を考える

そこで、配送ルート集合をR_0とした元の最適化問題の線形緩和問題を考えます。


\begin{aligned}
\mathrm{minimize}&\quad \sum_{r\in R_0}d_r x_r\\
\mathrm{subject\ to}&\quad \sum_{r\in R_0}a_{ri} x_r \geq 1,\quad i\in \mathrm{S}\setminus\{0\}, \\
&\quad 0\leq x_{r}\leq 1, \quad r\in R_0.\\
\end{aligned}

元の問題とその線形緩和問題の差分は、x_rの取る値の範囲だけですが、計算の効率は線形緩和問題の方が良いです。列生成法のアイディアは、「線形緩和した問題に解が改善しなくなるまで配送ルートを追加し、得られた配送ルート集合は、元の最適化問題の配送ルート集合としても、良いであろう」です。

線形緩和問題の双対問題を考える

最適化の世界では、「主問題が難しければ双対問題を考えろ」という言葉があるそうです。 線形緩和問題の双対問題を考えると、


\begin{aligned}
\mathrm{maximize}&\quad  \sum_{i \in S\setminus\{0\}}y_i\\
\mathrm{subject\ to}&\quad \sum_{i\in S\setminus\{0\}}a_{ri} y_i\leq d_r,\quad r\in R,\\
&\quad y_i\geq 0,\quad i\in \mathrm{S}\setminus\{0\}.\\
\end{aligned}

となります。上式の1つ目の制約条件から、線形計画緩和問題の各列(=変数)が双対問題の各行(=制約条件)に対応していることが分かります。線形計画問題で良い列の組合せを求めることは、双対問題では、新たな行(=制約条件)を追加しても最適解が変わらないような行の組合せを求めることに対応します。

列生成法の手順

列生成法の具体的な手順を示します。

  1. 初期配送ルート集合R_0を作ります。本記事では、倉庫から各店舗へピストン輸送するようなシンプルな配送ルート集合としました。
  2. 線形緩和問題の双対問題にR_0の元を制約式とし、双対問題の解y^{\ast}_iを得ます。
  3. 得られた解y_i^{\ast}の元で、容量制約をみたし、\displaystyle d_r^{\ast}-\sum_i a_{r^{\ast}i } y_i^{\ast}が負になるルートr^{\ast}を探します。
  4. ルートr^{\ast}を新たな制約条件として、R_0に追加したものをR_1とします。
  5. 手順2-4の手続きをR_0R_1に読み替えて実施します。
  6. 2-5の手続きを\displaystyle d_r-\sum_s a_{ri } y_iが負になるルートが見つからなくなるまで繰り返します。
  7. 得られた配送ルート集合R_kを元の最適化問題の配送ルート集合として考えて、解を求めます。

\begin{aligned}
\mathrm{minimize}&\quad \sum_{r\in R_k}d_r x_r\\
\mathrm{subject\ to}&\quad \sum_{r\in R_k}a_{ri} x_r \geq 1,\quad i\in \mathrm{S}\setminus\{0\},\\
&\quad x_{r}\in\{0, 1\},\quad r\in R_k.\\
\end{aligned}

手順6で、\displaystyle d_r^{\ast}-\sum_i a_{r^{\ast}i } y_i^{\ast}が負になるルートがないことは、双対問題視点で考えると、これ以上制約条件を追加しても、最適解が変化しないことを意味します。これは、主問題(=元の線形緩和問題)視点では、これ以上解の候補を追加しても、最適解が良くならないということです。つまり、線形緩和問題において、解の候補の列挙の完了を意味します。繰り返しになってしまいますが、列生成法は「線形緩和した問題に解が改善しなくなるまで解の候補を追加し、得られた解の候補集合は、元の最適化問題の解の候補集合としても、良いであろう」と考えます。

双対問題の制約条件を破る配送ルートの生成方法

列生成法の手順3についてもう少し説明致します。手順3では、


\begin{aligned}
d_r^{\ast}-\sum_{i \in S} a_{r^{\ast}i } y_i < 0
\end{aligned}

となる配送ルートr^{\ast}を見つける問題を考えます。この問題を列生成法の子問題と呼びます。子問題は、各店舗iに双対問題の解として得られたy_i^{\ast}で定義される利得がある前提で、倉庫から出発し、移動距離に対するペナルティを考慮しながら、各店舗を回り利得を集め、倉庫に戻る巡回路を求める問題とみなせます。第1項は、ルートr^{\ast}でかかる移動のコスト、第2項は、ルートr^{\ast}で訪れた店舗iから得られる利得を表します。この問題を、試しにMIPで定式化してみます。まず、店舗iから店舗jへ移動するかどうかを表す変数z_{ij}を用意します。


\begin{aligned}
z_{ij} = \begin{cases}1 \quad 店舗iから店舗jへ移動する\\0 \quad 店舗iから店舗jへ移動しない\end{cases}
\end{aligned}

この変数を用いた定式化は、


\begin{aligned}
\mathrm{minimize}&\quad \sum_{i!=j} c_{ij}z_{ij}-\sum_{i!=j} y_j^{\ast} z_{ij}\\
\mathrm{subject\ to}&\quad \sum_{i=1}^{N}z_{0i}=1,\\
&\quad \sum_{i=1}^{N}z_{i0}=1,\\
&\quad \sum_{j=1}^{N}z_{ij}\leq 1, \quad i=1,...,N,\\
&\quad \sum_{i=1}^{N}z_{ij}\leq 1, \quad i=1,...,N,\\
&\quad \sum_{j!=i}z_{ij}-\sum_{j!=i}z_{ji} =  0, \quad i=1,...,N,\\
&u_i+1-\mathrm{BigM}(1-z_{ij})\leq u_j, \quad i,j=1,...,N,\\
&\quad \sum_{j!=i}D_j z_{ij}\leq \mathrm{C},\\
&\quad 0\leq u_i \leq N,\quad u_i \in \mathbb{Z},\quad i=1,...,N,\\
&\quad z_{ij}\in\{0, 1\},\quad i!=j.\\
\end{aligned}

ただし、

  • C:トラックの積載容量
  • D_i:店舗iへ配送する荷物量
  • c_{ij}:店舗iと店舗jの移動距離
  • \mathbb{Z}:整数の集合
  • \mathrm{BigM}:十分大きな定数

制約式の意味ですが、

  • (1式目)倉庫から必ずトラックは1台出発するものとする。
  • (2式目)倉庫にどこかから必ずトラックが1台返ってくることとする。
  • (3式目)店舗iから出発するトラックは、最大でも1台とする。
  • (4式目)店舗jに到着するトラックは、最大でも1台とする。
  • (5式目)もし店舗iから出発するトラックがあるならば、店舗iへどこかの店舗から1台トラックが来ていることとする。*5
  • (6式目)倉庫以外の店舗から出発し、店舗を訪問し、倉庫以外の店舗に戻るようなトラックはないものとする。*6
  • (7式目)トラックは容量制約を満たすものとする。

となります。制約の6式目がわかりにくいかもしれません。考え方は、訪問する店舗iに訪問する順番で昇順にu_iという値を持つように制約条件を与えると、部分巡回路(=倉庫を含まない巡回路)が解として許されなくなるというものです。部分巡回路の取り扱い方について、参考文献[2]に詳しい記載がありますので良ければ、ご参照頂けたらと思います。子問題の解法は、上記のMIPの定式化に限らず、既存の解で制約条件を破る配送パターンを見つければ良いので、方法は様々あると思います。今回は、計算コストを考慮せず、定式化と実装の簡便さ、拡張性を意識して、MIPにて定式化してみました。

数値実験結果

サンプルデータについて、これまでに説明した定式化を解いてみます。下記のような条件で倉庫と店舗を用意します。

  • 配送店舗数:20
  • 倉庫の場所:(0.5,0.5)*7
  • 店舗の場所:(0,0),(0,1),(1,0),(1,1)の4点で囲まれる領域内で乱数
  • 店舗の物量:[5,15] *8を満たす整数の中で一様乱数
  • トラックの積載容量:100

pythonで実装し、最適化問題はPuLPを用いて解きました。
列生成法で子問題を解いた回数をステップ数と呼ぶことにします。ステップ数毎の子問題の目的関数の値は以下のようになりました。ステップ数が増えるにつれて、目的関数値が振動しながら0に近づいています。

f:id:bp-writer:20201020191247p:plain
ステップ数と子問題の目的関数値
ステップ数が159回目で、目的関数値が正になりました。つまり、158個の列生成により得られた解の候補の組み合わせで、列生成法による解を得ました。今回はトラックの台数が2台で、ステップ数=121,153回目に得られた解の候補が最終的に選択されました。

得られた配送計画は以下の通りです。

f:id:bp-writer:20201020190726p:plain
最適化問題の解として得られた配送計画

まとめ

配送計画問題を題材に、列生成法の説明をしてみました。わかりにくい点もあったかと思います。もし「自分でも列生成法を試してみたい」と思っていただけましたら幸いです。
最後にお約束ですが...プレインパッドは今、数理最適化に真剣に取り組んでいます。そして、こういった取り組みに興味のあるエンジニア・データサイエンティストを募集しています。数理最適化に興味のある方は、直近の転職意向の有無に関わらず気軽にお声がけください!

参考文献

  1. 宮本 裕一郎, はじめての列生成法, オペレーションズ・リサーチ, 57 (2012), 198-204.
  2. 久保 幹雄, ジョア・ペドロ・ペドロソ, 村松 正和 , アブドル・レイス , あたらしい数理最適化: Python言語とGurobiで解く, 近代科学社, 2012.

*1:半月前の私のような方を想定しています。

*2: 第1回のブログの題材がまさに定式化のTIPでした。第1回:最近学んだ数理最適化の定式化のチップスたち 【ブレインパッドの数理最適化ブログ】 - Platinum Data Blog by BrainPad

*3:参考文献[2]が非常に参考になると思います。

*4:本記事の列生成法の説明は参考文献[1]を参考にしています。筆者は半月前に初めて読んだ「はじめての列生成法」がわかりやすく、非常に面白かったので、文献とは別の最適化問題に対しても列生成法を試してみようと思い、本記事を執筆しました。

*5:この制約により巡回路を生成します。

*6:部分巡回路(=倉庫を含まない巡回路)を除去する制約でMTZ制約と呼ばれるものになります。

*7:ここで、(x,y)は、x,y座標を表します。

*8:ここで、[a,b]は閉区間を表します。