フレクトのクラウドblog re:newal

http://blog.flect.co.jp/cloud/からさらに引っ越しています

トップ > オペレーションズ・リサーチ > 港の在庫量を考慮した船舶の配送計画問題(Inventory Routing)を解く

港の在庫量を考慮した船舶の配送計画問題(Inventory Routing)を解く

オペレーションズ・リサーチ

こんにちは，研究開発室の福井です．

現在研究開発室においてオペレーションズ・リサーチ(OR)のビジネス活用について研究を行っております．OR の有名な問題として配送計画問題というものが存在しますが，海運における派生問題に Inventory Routing Problem (IRP)と呼ばれるものがあります．IRPでは，各港の在庫量に関する制約条件を守りつつ最適な船舶の航行計画を数理最適化により計画の出力を行います．本稿では，この IRP に関して解説を行いその後実際にモデリングをして例題を解いていきます．

はじめに: 海運における配送計画問題

海運ではその設定に応じて様々な計画問題を解く必要がありまが，ここでは Christiansen et al. 2013 における分類にもとづいて説明します．まず，貨物船の運用方法が定期船(liner ship)か不定期船(tramp ship)かによって分類することができます．定期船は決められた航路上を決められたスケジュールに従って航行し，その道中において貨物に荷積みと荷降ろしを行います．定期船では主にコンテナ等の輸送が行われており，コンテナ船は基本的には定期船に該当すると言えます．一方，不定期船では荷主の依頼等に従って貨物の輸送が行われ，粉や液体などのバルク貨物等が扱われることが多いようです．

次に定期船と不定期船のそれぞれの運用において必要となる計画問題について説明します．定期船で必要となる計画問題には最適な定期航路の設計と航路に対する船舶の割当を行う network design の問題(Agarwal & Ergun 2008)や fleet deployment (Gelareh & Meng 2010)等が存在します．一方，不定期船においてはコストを最小化(もしくは収益を最大化)するように，その時々で輸送しなければならない貨物に応じた各船舶の航路の決定と各船舶に関する荷積み/荷降ろしの計画作業(どの船舶に何をどれだけの量を積載してどこからどこへ運ぶか?)が必要となります．この不定期船の計画問題は設定に応じて更に 2 つに分類すことができ 1 つは cargo routing，もう 1 つは inventory routing と呼ばれます．Cargo routing では荷積み/荷降ろしの期間(time window)と貨物の詳細(量や貨物の内容など)が荷主側との調整で事前に決まっており，それらの指定された期間に港に訪問して荷積みと荷降ろしができるような配送計画を立てます*1．もうひとつの問題設定である inventory routing では荷積み/荷降ろしの期間や量が明示されておらずその代わりに各港が保有する製品や材料等の在庫量の推移(生産量/消費量)に関する情報が与えられます．Inventory routing のゴールは各港における在庫(inventory)を切らさない/溢れさせないよう船舶を用いて港間の在庫の輸送時期と輸送量を同時に計画する問題となります(図 1)．

本稿では，不定期船の計画問題のひとつである inventory routing についてその定式化から実際に計算機を用いて配送計画を出力し結果を確認してゆきます．

f:id:kazukifukui:20210608154806p:plain
図1 Inventory routingの概要．Inventory routingでは各港の在庫(赤いバー)を切らさない/溢れさせないように在庫の輸送を行います．

Inventory Routing Problem

ここでは改めて Inventory Routing Problem (IRP)の詳細について解説していきます．まず，IRP で以下のような情報が与えられます:

  • 港に関する情報の例
    • 各港間の距離
    • 計画開始時点での各港の初期在庫量
    • 各港での在庫の生産/消費レートの計画
    • 各港が保有できる在庫量の上限
    • 各港が保有する在庫量の安全な下限
  • 船舶に関する情報の例
    • 各船舶の航行速度
    • 各船舶の(距離などに応じた)航行コスト
    • 計画開始時点に各船舶に積載されている在庫量
    • 各船舶に積載できる在庫量の上限

これらの情報にもあるように，港の在庫はある量からスタートしその港において事前に定められている生産/消費レートに応じて港が持つ在庫を増加/減少させていくこととなります．当然ではありますが何もしなければ生産を行っている港では在庫が溢れてしまい，在庫を消費する港では在庫が底をついてしまいます．IRP はある生産拠点で材料を生産し別の拠点で加工や販売を行うような業務においてよく発生する問題設定となり，各拠点での在庫量が溢れたり底をついたりした場合業務に悪影響が出てしまいます．そのため IRP では各港での在庫量に関する制約条件を守りつつ，その上で最適な港の訪問順序(待ちを含む)と各港での在庫の荷積み/荷降ろし量の計画を行います(図 2)．

f:id:kazukifukui:20210608154850p:plain
図2 在庫量の消費と補充の例．

Time-Space Network による定式化

次に，IRP を数理計画問題として表現する際の定式化の方法について説明します．特に今回は IRP の定式化の 1 つである time-space network (Song & Furman 2013)を用いた定式化について説明します．Cargo routing (もしくは通常の車両の配送計画問題)では訪問先の港がネットワーク上の node で表され，node 間の edge に移動に要するコストや移動に必要な時間が設定されます．一方 IRP では各港において時間ごとに異なる在庫の生産/消費レートや時間毎に異なる在庫量を扱う必要が出てくるためある港を単一の node で表現するには不十分となります．そこで，time-space network では図 3 にあるように 1 つの港に対して時刻別に別々の node を用意します．図 3 における時間方向の node の個数は計画期間(Planning horizon)を表しており，例えば 1 ヶ月(30 日)分の計画を 1 日単位で離散的に扱う場合だと計画期間の長さは 30 というになります．Time-space network ではこの離散的に表された時刻のことを time period と呼びます．また，図における source node と sink node は計画の開始と終了を表しています．車両の配送計画の場合，車両はデポを出発してデポに返ってくるという設定が多いですが船舶においてはデポが存在せず海上や港などの任意の場所で計画が開始し任意の場所で計画が終了します．そのため，図中の source node と sink node は実際の港に対応しない仮想的な node となります．港間の移動時間は node 間に設定された arc (有向辺)で表されており，例えば港 i を時刻 t に出発して港 j への移動するさいの時間が(time period 換算で) 2 の場合，time-space network 上に node (i, t) から node (j, t + 2) への arc が追加されます．また，node (i, t) から (i, t + 1) への arc は港 i での船舶の待ちを表しています．

f:id:kazukifukui:20210608154920p:plain
図3 Time-space network．(Song & Furman 2013より引用)

次に，time-space network を用いた IRP の具体的な定式化を見ていきます．まずは，必要な記号を以下のように定義します．

  • 集合と index
    • t \in \mathcal{T} : Time period の集合
    • v \in \mathcal{V} : 船舶の集合
    • j \in \mathcal{J}^P : 在庫を生産している港の集合
    • j \in \mathcal{J}^C : 在庫を消費する港の集合
    • \mathcal{J} = \mathcal{J}^P \cup \mathcal{J}^C : 全ての港の集合
    • \mathcal{N} = \{ n=(j, t)\, |\, j \in \mathcal{J}, t \in \mathcal{T} \} : Time-space network 上の港に対応した node (港と time period のペア)の集合
    • \mathcal{N}_{s,t} : \mathcal{N}に source node n_{s} と sink node n_{t}を加えた全ての node の集合
    • a \in \mathcal{A} : Time space 上の全 arc の集合( \mathcal{A} \subset \mathcal{J} \times \mathcal{J})
    • a \in \mathcal{A}^v : 船舶 v に関連した arc の集合(船舶 v が通れる arc の集合)
    • a \in \mathcal{FS}^v_n : 船舶 v \in \mathcal{V} に関して node n \in \mathcal{N}_{s,t} から出ていく arc の集合
    • a \in \mathcal{RS}^v_n : 船舶 v \in \mathcal{V} に関して node n \in \mathcal{N}_{s,t} へと入っていく arc の集合
  • パラメータ
    • B_j : 港 j \in \mathcal{J} でのバースの容量(停泊可能な船舶の数)
    • C^v_a : 船舶 v が arc a \in \mathcal{A}^v を航行した際に発生するコスト
    • d_{j,t} : 港 j \in \mathcal{J}，time period t \in \mathcal{T} における在庫の生産/消費量
    • \Delta_j : 港における在庫の生産/消費を表す符号． j \in \mathcal{J}^P なら +1， j \in \mathcal{J}^Cなら -1
    • Q^v : 船舶 v の容量
    • S^{max}_j : 港 j \in \mathcal{J} の容量
    • s_{j,0} : 港 j \in \mathcal{J} の在庫量の初期値
    • s^v_0 : 船舶 v \in \mathcal{V} の在庫量の初期値
  • 決定変数(最適化によって決まる値)
    • f^v_n : 船舶 v \in \mathcal{V} が node n \in \mathcal{N} にて荷積み/荷降ろしする在庫の量
    • s_{j,t} : Time period t \in \mathcal{T} の開始時点において港 j \in \mathcal{J} が保有する在庫量
    • s^v_t : Time period t \in \mathcal{T} の開始時点において船舶 v \in \mathcal{V} が保有する在庫量
    • x^v_a : 船舶 v \in \mathcal{V} が arc a \in \mathcal{A}^vを通るなら 1，そうでないなら 0
    • z^v_n : 船舶 v \in \mathcal{V} が node n \in \mathcal{J} において荷積み/荷降ろしを行うなら 1，そうでないなら 0

次に，上記の記号を用いて IRP の定式化を以下のように行います:

\displaystyle \begin{align} \min &\sum_{v \in \mathcal{V}} \sum_{a \in \mathcal{A}^v} C^v_a x^v_a \tag{1} \\ \text{s.t.} & \sum_{a \in \mathcal{FS}^v_{n_s}} = 1, & \forall v \in \mathcal{V} \tag{2} \\ & \sum_{a \in \mathcal{RS}^v_{n_t}} = 1, & \forall v \in \mathcal{V} \tag{3} \\ & \sum_{a \in \mathcal{FS}^v_{n}} x^v_a = \sum_{a \in \mathcal{RS}^v_{n}} x^v_a, & \forall v \in \mathcal{V}, \forall n \in \mathcal{N} \tag{4} \\ & s_{j,t} = s_{j,t-1} + \Delta_j \left(d_{j,t-1} - \sum_{v \in \mathcal{V}} f^v_{(j, t-1)} \right), & \forall n = (j,t) \in \mathcal{N} \tag{5} \\ & s^v_t = s^v_{t-1} + \sum_{j \in \mathcal{J}} \Delta_j\, f^v_{(j, t-1)}, & \forall t \in \mathcal{T}, \forall v \in \mathcal{V} \tag{6} \\ & \sum_{v \in \mathcal{V}} z^v_n \le B_j, & \forall n = (j,t) \in \mathcal{N} \tag{7} \\ & z^v_n \le \sum_{a \in \mathcal{RS}^v_n} x^v_a, & \forall n \in \mathcal{N}, \forall v \in \mathcal{V} \tag{8} \\ & 0 \le s_{j,t} \le S^{max}_j & \forall j \in \mathcal{J}, \forall t \in \mathcal{T} \tag{9} \\ & 0 \le s^v_t \le Q^v, & \forall v \in \mathcal{V}, \forall t \in \mathcal{T} \tag{10} \\ & 0 \le f^v_n \le \min\{Q^v, S^{max}_j\} z^v_n, & \forall n=(j,t) \in \mathcal{N}, \forall v \in \mathcal{V} \tag{11} \\ & x^v_a \in \{0, 1\}, & \forall v \in \mathcal{V}, \forall{a} \in \mathcal{A}^v \tag{12} \\ & z^v_n \in \{0, 1\}, & \forall n \in \mathcal{N}, \forall v \in \mathcal{V} \tag{13} \end{align}


これらの数式の意味は以下のようになります．

  1. 数理計画問題の目的関数．全船舶に関する航行コストの合計を最小化する
  2. 船舶は必ず(仮想の) source node から出発する
  3. 船舶は計画終了時に必ず(仮想の) sink node へ訪れる
  4. Node 上の移動に関する整合性．船舶が node に入ったら必ずその node から出ていく
  5. 港の在庫量の推移に関する整合性
  6. 船舶に積載された在庫量の推移に関する整合性
  7. 港にて同時に荷積み/荷降ろしする船舶数がその港のバースの数を超えない
  8. 訪問した港でのみ荷積み/荷降ろしが可能
  9. 港の在庫量が上限/下限の範囲に収まっている
  10. 船舶の在庫量が上限/下限の範囲に収まっている
  11. 荷積み/荷降ろしする際のみ荷積み/荷降ろし量が正の値を取る

ここで紹介した定式化は必要最低限の IRP の設定であり，実際には業務内容やルールに応じて以下のような制約条件が追加で考慮される場合があります．

  • 在庫の荷降ろしによって発生する収益の考慮
  • 港での在庫の保管に伴うコストの考慮
  • 喫水制限の考慮
  • メンテナンスのためのドック入りの考慮
  • 運搬中の在庫のボイルオフ(LNG，LP ガスなどの気化)の考慮
  • 配送計画と生産量の計画との同時最適化

例題: MIRPLIB

ここでは実際に IRP の問題を実装し OSS の汎用ソルバを用いて計画問題を解いた結果をお見せします．先にも述べたように業務に応じて様々な制約条件が発生しますがここでは IRP のベンチマークの 1 つである MIRPLIB (Papageorgiou+ 2014)を使用します．

MIRPLIB では規模や設定が異なる IRP のインスタンス(問題例)を複数提供しており，以降ではその内の 1 つを例として用いていきます*2．図 4 は例題で使用するインスタンスで扱う港のx-y座標を示したものです．図に示されているように， 1 ヶ所の loading port (在庫を生産しており船舶が荷積みを行う港)と 3 ヶ所の discharging port (在庫を消費しており船舶が荷降ろしを行う港)が存在します．また，このインスタンスでは 7 隻の貨物船(速度，容量は全て同じ)が利用可能となっています．

図4 例題の港の座標と位置関係．(カーソルを重ねると詳細が見られます．)

OR-Tools を用いた IRP の実装

問題の実装にはOR-Tools の CP-SAT ソルバを使用します．OR-Tools は Google が開発した OR のための OSS のツールであり，様々な問題に対するソルバを提供しています．以下にOR-Toolsを用いた簡単な制約付き最適化問題の実装例(Python)を記載します(下記はOR-Toolsのプロジェクトページ上の公式のサンプルとほぼ同様のものとなります)．

from ortools.linear_solver import pywraplp

# BackendにCP-SATを指定してソルバを用意
solver = pywraplp.Solver.CreateSolver('CP-SAT')

# 非負の整数の変数x, yを用意
infinity = solver.infinity()
x = solver.IntVar(0.0, infinity, 'x')
y = solver.IntVar(0.0, infinity, 'y')

# 制約条件を設定
# x + 7 * y <= 17.5
solver.Add(x + 7 * y <= 17.5)
# x <= 3.5
solver.Add(x <= 3.5)

# 目的関数を設定(ここでは最大化)
solver.Maximize(x + 10 * y)

# ソルバによる解の探索を実行
status = solver.Solve()

上記の例ように制約条件と目的関数を指定することで様々な数理計画問題を解くことが可能となり，式 1 - 式 13 の制約条件と目的関数に関しても同様に 1 つずつプログラムに書き起こしていくことでIRPを実装することができます*3 *4．

結果

OR-Tools を用いて今回の例題を解くことで得られた各港と船舶上の在庫量の時間変化が図 5，6 となります．得られた計画において，各 discharging port の在庫量が尽きる前に船舶により在庫の補充が行われていることがわかります．また，図 7 は計画において各船舶が港を訪問する順序をしめしており，図 6 と図 7 より出力された計画において各港のバースの制限が守られていることがわかります．また，MIRPLIB の例題に少し変更を加えた際の解への変化を見るために制約条件を追加してみます．ここでは，何らかの理由により time period が12 ~ 19 の間は DischargeRegion_0_Port_2 が利用不可となる状況を想定しています．これを表現する制約条件は至ってシンプルで，time period が 12 ~ 19 に対応する node に入ってくる arc を通る船舶が 0 隻になるように設定することで実現できます．図 8，9 に示した結果を見ると，設定した通り特定の time period の間は船舶が DischargeRegion_0_Port_2 に訪れていないことがわかります．また，その期間内は在庫の補充が行えないため，アクセス不可となる期間の開始直前と直後に在庫が補充されるような計画が出力されていることがわかります．

下図では(あくまで私個人の感想ですが)複雑な計画が出力されており，人間が実際に同様の計画作業を行う場合は恐らく長い時間を要するものと考えられます．しかし，配送計画の作業を OR の問題として計算機に解かせることで人手に頼らずに計画の候補*5を得ることができました*6．

図5 各港の在庫量．(カーソルを重ねると詳細が見られます．また，港名をクリックすることで表示の on/off が切り替えられます．)
図6 各船舶の積載量．(カーソルを重ねると詳細が見られます．また，船舶名をクリックすることで表示の on/off が切り替えられます．)
図7 航行スケジュール．(カーソルを重ねると詳細が見られます．また，船舶名をクリックすることで表示の on/off が切り替えられます．)
図8 各港の在庫量．時刻 [12, 19] にて DischargeRegion_0_Port_2 へ訪問不可とした場合．
図9 航行スケジュール．時刻 [12, 19] にて DischargeRegion_0_Port_2 へ訪問不可とした場合．

まとめ

本稿では船舶の配送計画問題の 1 つである inventory routing problem について解説し，実際に MIRPLIB のインスタンスの 1 つを例題として解いてその結果を示しました．例題の結果より，人間ならば時間を要するような計画作業でも計算機上で汎用ソルバを用いることで計画の候補が得れることがわかりました．今回の例題ほどの規模なら比較的早く結果が得られますが，開発室内の実験では港数等の増加に伴い計算に要する時間が大幅に増加する事がわかっており，規模の大きな問題に対する高速化が今後の課題となります．また，現実においては discharging port 側での在庫の消費量が市場での需要量に応じて変化する場合があり，その際には discharging port での消費量が確率的な不確かさを持つためそのような不確かさの考慮したモデリングも今後課題となります．

最後に

最後まで読んでいただきありがとうございました．フレクトでは今後もお客様にとって付加価値となるような新たな分野や技術に関して開拓を行っていきます．今回紹介した船舶の配送計画問題以外の様々な OR の諸問題に関しても研究を行っているため，もしご興味がありましたら是非ご相談ください．

参考文献

  • Agarwal, R., & Ergun, Ö. (2008). Ship scheduling and network design for cargo routing in liner shipping. Transportation Science, 42(2), 175-196.
  • Christiansen, M., Fagerholt, K., Nygreen, B., & Ronen, D. (2013). Ship routing and scheduling in the new millennium. European Journal of Operational Research, 228(3), 467-483.
  • Gelareh, S., & Meng, Q. (2010). A novel modeling approach for the fleet deployment problem within a short-term planning horizon. Transportation Research Part E: Logistics and Transportation Review, 46(1), 76-89.
  • Papageorgiou, D. J., Nemhauser, G. L., Sokol, J., Cheon, M. S., & Keha, A. B. (2014). MIRPLib–A library of maritime inventory routing problem instances: Survey, core model, and benchmark results. European Journal of Operational Research, 235(2), 350-366.
  • Song, J. H., & Furman, K. C. (2013). A maritime inventory routing problem: Practical approach. Computers & Operations Research, 40(3), 657-665.

*1:Cargo routing はトラックなどによる陸上での(pickup & deliveryの)配送計画問題とほぼ同じように定式化することができます

*2:具体的には Group 1 の LR1_1_DR1_3_VC1_V7a を planning horizon=30 で区切って使用します

*3:IRP は本来混合整数計画問題ですが，今回使用する CP-SAT の仕様の関係で整数計画問題としてモデリングしています．その際，小数点以下の値が発生する数式に関しては全体に大きな整数を掛けることで小数点以下の数字を除去しています．

*4:厳密には，MIRPLIB のインスタンスでは式 1-式 13 に加えてもういくつかの制約条件が追加されます．詳細に関しては Papageorgiou+ 2014 をご確認ください

*5:ここで"候補"と記載しているのは，得られた解が最適かどうか証明できていないからです．CP-SAT では全ての解を探索し終わるかもしくは現在の解が解の bound と一致した際にその解が最適であると判定しており，それ以外の状況で得られた解に関しては「最適かどうかはわからない実行可能解の1つ」として出力されます．今回の例題で得られた解は後者に該当するためあくまでも実行可能な計画の候補の1つということになりますが，もし最適性を証明したい場合はより長い時間の間計算を回す必要があります．

*6:今回の例題に関してなんらかの実行可能解を得るだけなら Core i5 の MacBook Pro を用いて 5 分程で済みます