組合せ最適化でコンテナ納品計画を作成してみた

#はじめに

斉藤努先生の記事[「組合せ最適化でナーススケジューリング問題を解く」](https://qiita.com/SaitoTsutomu/items/a33aba1a95828eb6bd3f)を以下で読み替えて、倉庫へのコンテナ納品計画最適化に応用してみました。

個人的には（詳細を詰めれば）実務でも利用できそうだと感じたので備忘録も兼ねて残します。
*内容はありふれたものですので、中～上級者の方は読み飛ばしてください。

|#  |読替え前  |読替え後  |
|---|---------|---------|
|1  |従業員    |コンテナ  |
|2  |必要人数  |倉庫受入枠  |

#前提条件
* 1つの会社で複数部署が各商材をコンテナで輸入している。
* 日本の港にコンテナが到着した後、近隣の倉庫へ納品する。

- 倉庫納品は事前に予約が必要（2週目と4週目に翌2週間分を予約）。予約方法は以下の通り。
 - 各部署の担当者が自部署の都合の良い納品日（納品希望日）をピンポイントで指定して倉庫へ連絡する。
 - 倉庫で受入枠（納品可能なコンテナ本数）は平日：最大5本/日、土日：0本/日。
 - 倉庫は各部署の希望日と受入枠を勘案して、納品計画を作成する。
  - 各希望日で倉庫受入枠≧合計本数の時は、全数希望日で予約が確定する。
  - 各希望日で倉庫受入枠＜合計本数の時は、超過分は別日で予約が確定する。
  - 希望日以外に予約を確定せざるを得ない場合は、倉庫は都度各部署と調整する必要がある。


#考えたこと
- ピンポイントで納品希望日に納品しなくても、その日までに倉庫に商材があればよいのではないか。
- すなわち、入港日翌日～納品希望日までのいずれかの日に納品できれば問題ないのではないか。（以下イメージ）
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/985685/57e28a1a-6dda-d9e4-adc0-bafc4fd96d18.png)



#やりたいこと
 - pythonを用いて、なるべくすべてのコンテナを入港日翌日～納品希望日の間に納品するような組合せを求めたい。
 - 希望日より納品が早まっても困ることはない（とする）。むしろなるべく早く倉庫へ納品する方がよい。
 - 一方で、客先への納期もあるので、希望日から遅れて納品するのはなるべく避けたい。
 -  ただし、入港日以前まで納入を早めることは物理的に不可。

#pulpの利用
- 引用記事同様に無料で使える最適化ライブラリとしてpulpを使用。
- pulpの補助ライブラリとしてortoolpyも使用。

```terminal:terminal
$ pip install pulp
$ pip install ortoolpy
```

#コンテナ予約表のサンプル
- 各部署からの納品希望日を以下devan_booking.csvとして一元化している前提。
- ランダムで起こした以下の2/1～2/15納品希望分でシミュレーション。

|コンテナ番号|部署|入港日|希望納品日|
|:----|:----|:----|:----|
|コンテナA01|A部署|2022/2/5|2022/2/11|
|コンテナA02|A部署|2022/2/3|2022/2/6|
|コンテナA03|A部署|2022/2/3|2022/2/10|
|コンテナA04|A部署|2022/1/31|2022/2/3|
|コンテナA05|A部署|2022/1/31|2022/2/7|
|コンテナB01|B部署|2022/2/9|2022/2/15|
|コンテナB02|B部署|2022/1/31|2022/2/2|
|コンテナB03|B部署|2022/2/4|2022/2/7|
|コンテナB04|B部署|2022/2/5|2022/2/11|
|コンテナB05|B部署|2022/1/31|2022/2/2|
|コンテナC01|C部署|2022/1/27|2022/2/1|
|コンテナC02|C部署|2022/2/2|2022/2/8|
|コンテナC03|C部署|2022/2/3|2022/2/10|
|コンテナC04|C部署|2022/1/31|2022/2/2|
|コンテナC05|C部署|2022/2/8|2022/2/14|
|コンテナD01|D部署|2022/1/26|2022/2/2|
|コンテナD02|D部署|2022/2/3|2022/2/5|
|コンテナD03|D部署|2022/2/5|2022/2/8|
|コンテナD04|D部署|2022/2/5|2022/2/11|
|コンテナD05|D部署|2022/1/26|2022/2/2|
|コンテナE01|E部署|2022/2/9|2022/2/12|
|コンテナE02|E部署|2022/2/1|2022/2/7|
|コンテナE03|E部署|2022/2/2|2022/2/8|
|コンテナE04|E部署|2022/2/10|2022/2/15|
|コンテナE05|E部署|2022/2/9|2022/2/13|
|コンテナF01|F部署|2022/2/9|2022/2/11|
|コンテナF02|F部署|2022/1/27|2022/2/3|
|コンテナF03|F部署|2022/2/4|2022/2/10|
|コンテナF04|F部署|2022/1/26|2022/2/1|
|コンテナF05|F部署|2022/2/9|2022/2/11|
|コンテナG01|G部署|2022/2/4|2022/2/9|
|コンテナG02|G部署|2022/2/1|2022/2/8|
|コンテナG03|G部署|2022/1/30|2022/2/6|
|コンテナG04|G部署|2022/2/5|2022/2/11|
|コンテナG05|G部署|2022/2/7|2022/2/13|
|コンテナH01|H部署|2022/2/5|2022/2/8|
|コンテナH02|H部署|2022/1/31|2022/2/5|
|コンテナH03|H部署|2022/1/29|2022/2/5|
|コンテナH04|H部署|2022/1/29|2022/2/2|
|コンテナH05|H部署|2022/2/3|2022/2/10|
|コンテナI01|I部署|2022/2/6|2022/2/8|
|コンテナI02|I部署|2022/2/5|2022/2/8|
|コンテナI03|I部署|2022/2/3|2022/2/8|
|コンテナI04|I部署|2022/2/8|2022/2/10|
|コンテナI05|I部署|2022/2/6|2022/2/8|
|コンテナJ01|J部署|2022/1/29|2022/2/5|
|コンテナJ02|J部署|2022/2/5|2022/2/8|
|コンテナJ03|J部署|2022/2/6|2022/2/11|
|コンテナJ04|J部署|2022/2/9|2022/2/12|
|コンテナJ05|J部署|2022/1/31|2022/2/4|

#pythonで作ってみた
- 書き慣れていないので汚いコードですがご容赦頂ければ幸いです。

###ライブラリのインストール
```python:python
import numpy as np
import pandas as pd
from datetime import datetime as dt
from datetime import timedelta as td
from pulp import*
from ortoolpy import addvars, addbinvars
```

###コンテナ納品予約表(CSV)の読み込み
- CSVを読み込み、コンテナ数×日付の表へ変形。
- この後の計算の都合上、tbl2のTrue/Falseを入れ替え(tbl_2)。

```python:python
tbl1 = pd.read_csv('devan_booking.csv') #csvの読み込み

a = pd.to_datetime(tbl1['入港日'])
b = pd.to_datetime(tbl1['希望納品日'])
daymin = a.min()+ td(days=1) #すべてのコンテナの入港日のうち一番若い日付+1日を取得
daymax = b.max()  #すべてのコンテナの希望納品日のうち一番先の日付を取得
day = (daymax - daymin).days + 1 #日数
daylist =[daymin.date() + td(days = i) for i in range(day)] 
# 今回のサンプルはdaymin=1/27、daymax=2/15なので1/27-2/15の日付リストを作成する。
ctn = len(tbl1) #コンテナ数
tbl2 = np.zeros(ctn*day).reshape(ctn,day)
for i in range(ctn):
    wnt = (b[i]-daymin).days
    pod = (a[i]-daymin).days
    tbl2[i,pod+1:wnt]=1 # 入港日+1日～納品希望日まで1にする。
tbl2 = tbl2.astype(np.bool).copy() #希望日付を1に変更

tbl_2 = ~tbl2 #最小化問題のため、Trueで納品できたときは目的関数の計算で0になるようにする。~でFalseとTrueを裏返す。
tbl_2 = tbl_2.astype(np.int) #0と1に戻す。
```
###倉庫受入枠(tbl3)の作成
- 平日は5本/日、土日は0本/日とする。今回、祝日は考えないことにする。

```python:python
# 倉庫受入枠のテーブル。平日は5本/日、土日は0本/日とする。
tbl3 = np.zeros(day)
count = daymin
for i in range(day):
    if count.weekday() == 5 or count.weekday() == 6: #5が土曜、6が日曜
        tbl3[i] = 0
    else:
        tbl3[i] = 5
    count += td(days=1)
```
###入港前納品のペナルティ表(tbl4)の作成
- 入港前の納品は物理的に不可なので、そのような解が出ないよう特別なペナルティ(10000点)を課す。
- 入港後、速やかに納品された方が良いと考え、納品日が早い方が点数が良くなるようにした。

```python:python
tbl4 = tbl_2*10000
for i in range(ctn):
    pod = (a[i]-daymin).days
    count = 0 
    for j in range(day-pod):
        tbl4[i,j+pod]=count
        count+=1
    
# 入港日+1日以前には納品できないので、最適解にならないよう10000点を置く。
# 入港日+1日以降の納品は物理的には可能なので、0から数字を置いていく。
# なるべく前の日付で納品されるよう、若い順に0から点数をつける。
```

###変数の宣言
- コンテナ納品有無が入るコンテナ数×日付のバイナリ変数を作る(tbl5)。
- 受入枠より本数超過したにペナルティを与えるように納品本数-受入枠>0のときに数値が入る変数を作る。
 - 納品本数-受入枠≦0でも問題ないので＋の時のみ入るよう非負にする。

```python:python
tbl5 = np.array(addbinvars(day,ctn))
#コンテナ本数×日数の変数matrix方を作成(バイナリ変数)

tbl6 = addvars(day) # 納品オーバーを計算し入れるための変数（非負変数）を日数分作成

tbl_2 = tbl_2.T.tolist() #目的関数作成用に変換
tbl3 = tbl3.T.tolist()
tbl4 = tbl4.T.tolist()
```
###モデル化・最適化処理
- 条件を式に落とし最適化計算。
- ペナルティの点数は仮置き。今後改良の余地あり。

```python:python
Cwhs = 100 #受入枠を超えたときのペナルティ
Cnwnt = 10 #納品希望日以外のペナルティ(この定数とtbl4の内積がペナルティ)
Clt = 1 #リードタイムが1日増えると1点ペナルティ。LT重みは1のため不要だが、今後の修正のために残す。

m = LpProblem(sense=LpMinimize) #最小化問題の宣言
m += (Cwhs * lpSum(tbl6)
    + Cnwnt * lpDot(tbl_2,tbl5)
    + Clt * lpDot(tbl4,tbl5))
# 目的関数の式
# 1行目は納品制限本数を超過した本数当たり、納品オーバー分のペナルティ
# 2行目は納品可能日ではない日に納品した本数あたり、希望不可分のペナルティ
# 3行目は各コンテナ納品した日付でのLT合計。
# すべての合計が最小になるときのスケジュール表を出力

for i in range(day):
    m += tbl6[i] >= (lpSum(tbl5[i,:])-tbl3[i])
# 制約条件の式1
# 納品オーバーの計算式（納品本数ー最大納品数）になる。
# 納品本数オーバーは非負の変数で（納品本数ー最大納品数）以上の任意の正の値
# 目的関数上、任意の値のうち最小の値、すなわち（納品本数ー最大納品数）となる
# 2行目のーの式を入れたほうが納品本数のばらつきが収まる。

for j in range(ctn):
    m += lpSum(tbl5[:,j]) <= 1
    m += lpSum(tbl5[:,j]) >= 1

# 制約条件の式2
# 1コンテナあたりの納品日を必ず1つ決める

m.solve()
print('目的関数', value(m.objective))
```

###結果の出力(コード)
- 計算結果をCSVに出力。
- 最初に読み込んだdevan_booking.csvの右に最適化の結果を追記して、新たにdevan_opt.csvとして保存。

```python:python
result = np.vectorize(value)(tbl5).astype(int).T #ソルバーの出力結果（日にち×コンテナ数で0-1が埋まっている）
daykey = [i for i in range(day)] 
datedic = {key: date for key, date in zip(daykey, daylist)} 
result = [datedic[sum(result[i]*daykey)] for i in range(ctn)] #intからdateに変換

tbl1['入港日'] = a.copy()
tbl1['希望納品日'] = b.copy()
tbl1['最適納品日']= result
tbl1.to_csv('devan_opt.csv',index=False,encoding='utf-8_sig')   
```

|コンテナ番号|部署|入港日|希望納品日|最適納品日|希望範囲内か|
|:----|:----|:----|:----|:----|:----|
|コンテナA01|A部署|2022/2/5|2022/2/11|2022/2/8|TRUE|
|コンテナA02|A部署|2022/2/3|2022/2/6|2022/2/4|TRUE|
|コンテナA03|A部署|2022/2/3|2022/2/10|2022/2/4|TRUE|
|コンテナA04|A部署|2022/1/31|2022/2/3|2022/2/1|TRUE|
|コンテナA05|A部署|2022/1/31|2022/2/7|2022/2/2|TRUE|
|コンテナB01|B部署|2022/2/9|2022/2/15|2022/2/10|TRUE|
|コンテナB02|B部署|2022/1/31|2022/2/2|2022/2/1|TRUE|
|コンテナB03|B部署|2022/2/4|2022/2/7|2022/2/8|FALSE|
|コンテナB04|B部署|2022/2/5|2022/2/11|2022/2/9|TRUE|
|コンテナB05|B部署|2022/1/31|2022/2/2|2022/2/1|TRUE|
|コンテナC01|C部署|2022/1/27|2022/2/1|2022/1/28|TRUE|
|コンテナC02|C部署|2022/2/2|2022/2/8|2022/2/3|TRUE|
|コンテナC03|C部署|2022/2/3|2022/2/10|2022/2/8|TRUE|
|コンテナC04|C部署|2022/1/31|2022/2/2|2022/2/1|TRUE|
|コンテナC05|C部署|2022/2/8|2022/2/14|2022/2/11|TRUE|
|コンテナD01|D部署|2022/1/26|2022/2/2|2022/1/27|TRUE|
|コンテナD02|D部署|2022/2/3|2022/2/5|2022/2/4|TRUE|
|コンテナD03|D部署|2022/2/5|2022/2/8|2022/2/9|FALSE|
|コンテナD04|D部署|2022/2/5|2022/2/11|2022/2/8|TRUE|
|コンテナD05|D部署|2022/1/26|2022/2/2|2022/1/27|TRUE|
|コンテナE01|E部署|2022/2/9|2022/2/12|2022/2/10|TRUE|
|コンテナE02|E部署|2022/2/1|2022/2/7|2022/2/2|TRUE|
|コンテナE03|E部署|2022/2/2|2022/2/8|2022/2/3|TRUE|
|コンテナE04|E部署|2022/2/10|2022/2/15|2022/2/11|TRUE|
|コンテナE05|E部署|2022/2/9|2022/2/13|2022/2/11|TRUE|
|コンテナF01|F部署|2022/2/9|2022/2/11|2022/2/10|TRUE|
|コンテナF02|F部署|2022/1/27|2022/2/3|2022/1/28|TRUE|
|コンテナF03|F部署|2022/2/4|2022/2/10|2022/2/9|TRUE|
|コンテナF04|F部署|2022/1/26|2022/2/1|2022/1/27|TRUE|
|コンテナF05|F部署|2022/2/9|2022/2/11|2022/2/10|TRUE|
|コンテナG01|G部署|2022/2/4|2022/2/9|2022/2/8|TRUE|
|コンテナG02|G部署|2022/2/1|2022/2/8|2022/2/2|TRUE|
|コンテナG03|G部署|2022/1/30|2022/2/6|2022/1/31|TRUE|
|コンテナG04|G部署|2022/2/5|2022/2/11|2022/2/9|TRUE|
|コンテナG05|G部署|2022/2/7|2022/2/13|2022/2/11|TRUE|
|コンテナH01|H部署|2022/2/5|2022/2/8|2022/2/7|TRUE|
|コンテナH02|H部署|2022/1/31|2022/2/5|2022/2/2|TRUE|
|コンテナH03|H部署|2022/1/29|2022/2/5|2022/1/31|TRUE|
|コンテナH04|H部署|2022/1/29|2022/2/2|2022/1/31|TRUE|
|コンテナH05|H部署|2022/2/3|2022/2/10|2022/2/4|TRUE|
|コンテナI01|I部署|2022/2/6|2022/2/8|2022/2/7|TRUE|
|コンテナI02|I部署|2022/2/5|2022/2/8|2022/2/7|TRUE|
|コンテナI03|I部署|2022/2/3|2022/2/8|2022/2/4|TRUE|
|コンテナI04|I部署|2022/2/8|2022/2/10|2022/2/9|TRUE|
|コンテナI05|I部署|2022/2/6|2022/2/8|2022/2/7|TRUE|
|コンテナJ01|J部署|2022/1/29|2022/2/5|2022/1/31|TRUE|
|コンテナJ02|J部署|2022/2/5|2022/2/8|2022/2/7|TRUE|
|コンテナJ03|J部署|2022/2/6|2022/2/11|2022/2/10|TRUE|
|コンテナJ04|J部署|2022/2/9|2022/2/12|2022/2/11|TRUE|
|コンテナJ05|J部署|2022/1/31|2022/2/4|2022/2/1|TRUE|

- 2本だけ希望範囲内に収まなかった。（希望日から1日遅れ）

###結果の出力(グラフ)と効果の考察

- 上のCSVデータを棒グラフにすると以下のような結果となった。
 - 最適化処理前は何日か倉庫受入枠を超過しているが最適化処理後はすべて倉庫受入枠に収まっている。
 - 最適化処理前後で1日の納品本数が安定している。倉庫で人員確保がしやすい、手待時間が発生しにくい。
 - 最適化処理前後で入港～倉庫納品までのリードタイムが短くなっていそう。デマレージ・ディテンションチャージの発生抑制。客先からの急な需要増へ対応可、販売機会損失の抑制。C/F良化。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/985685/460fd4a7-26be-7a05-9b3d-c5a2c266c8d2.png)


#まとめ・感想
上記内容で業務にある程度までは利用できそうだと思いました。
具体的な金額は伏せますが、自社のデータを使い、同じシミュレーションをしたところ、ロスコスト抑制額、C/F良化額ともに驚くような金額がはじき出されました。
ただし、各部署の要望通り納品されない頻度が高くなりますので、関連組織と業務の再設計が必要かと思います。
ぱっと思いつく限り、以下の課題があります。

- 1/27～1/31までも通常の倉庫受入数(5本/日)で計算したが、実際にはこの日程は前回シミュレーションの範囲に含まれるので、既に予約確定しているコンテナがあるはず。5本は入らない。
- 各部署で月末在庫を管理している場合は、月初に納品希望を出したコンテナを前月末に納品することに反対される場合がある。
- 予約時点で最適化をしても、実際のコンテナ輸送時際に船が遅延して納品計画を一から組み直さないといけない場合がある(2022年1月時点では船の遅延率がかなり高いエリアもある)。
 
業務・運用面での課題はあるにしても、このpulpを使って誰でも無料で最適化計算ができるのは素晴らしいことだと思います。現場では依然、人海戦術で乗り切っている業務が多々ありますので、活用できる場は多々ありそうです。他に利用できる業務がないか探してみたいと思います。

#参考文献
- 引用記事の更に引用記事です。[「遺伝的アルゴリズムでナーススケジューリング問題（シフト最適化）を解く」](https://qiita.com/shouta-dev/items/1970c2746c3c30f6b39e)