GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(2)輸送最適化問題

#はじめに

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」の練習問題を解く第二回目の記事です。

第一回目の記事は以下です。まずはこちらをご覧下さい。
>GoogleのOR-Toolsで数理最適化モデルの練習問題を解く(1)一番易しいマス埋め問題
>https://qiita.com/ttlabo/private/7e8c93f387814f99931f

#輸送最適化問題

参考テキストに出てくる練習問題の2番目です。
以下の問題をさっそく行ってみましょう。

>:speech_balloon: 問題

```ruby:7.2.問題
倉庫群から工場群へ部品を搬送したい。総輸送費が最小となる計画を求めたい。各倉庫からの搬出は供給可能量以下とし、各工場への搬入は需要量以上とする。

■モデル
各倉庫から各工場への輸送量が、意思決定対象すなわち変数となる。目的関数は、変数に輸送費をかけて和を取った総輸送費となる。
守らなければならないことは、供給可能量を超えないことと需要量を満たすことである。
　変数　　　倉庫群から工場群への総輸送量
　目的関数　総輸送費→最小化
　制約条件　各倉庫からの搬出は、供給可能量以下
　　　　　　各工場への搬入は、需要量以上
```

>:question: **考え方**

ここではテキストの例に倣い、各倉庫(W1,W2,W3)の供給可能量と各工場(F1,F2,F3,F4)の需要量は以下とします。

|倉庫|供給可能量|
|:--:|:--:|
|W1|35|
|W2|41|
|W3|42|

|工場|需要量|
|:--:|:--:|
|F1|28|
|F2|29|
|F3|31|
|F4|25|

また、各倉庫から各工場へのそれぞれの(個別の)輸送費を以下の縦横表[[表1]](マトリックス表)の通りとします。部品1個あたりの輸送費です。
![001.JPG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/4d2c98bb-0e8b-6546-d1df-3e439e0d21a4.jpeg)

一方、同様に輸送量について下表[[表2]](マトリックス表)の通りとします。輸送する部品の個数です。
縦横表となっているため、それぞれの座標(行,列)における輸送量をamount(i, j)とします。
iは行位置を示し、1行～3行まであることからそれぞれの値を0,1,2とします。
jは列位置を示し、1列～4列まであることからそれぞれの値を0,1,2,3とします。

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/c5022ac4-c4df-9232-9d3c-bac2014af3ce.jpeg)


総輸送量はそれぞれの輸送費と輸送量を掛けて総和を取ったものとなります。
つまり、表1と表2のそれぞれのセルの値同士をかけてから総和を取ったものです。


総輸送量が最小となる数理モデルを考え、解を求めます。
すなわち、目的関数が最小となる最小化モデルの混合整数計画問題として考えます。



>:a: **解答**

練習問題の第一回目と同様におまじないを定義します。

```7.2.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp
```

混合整数計画問題としてCBCソルバーを使用することを宣言します。

```7.2.renshu.py
# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)
```

表1に従い、輸送費の定義を行います。
ここでは、2次元配列で定義します。

```7.2.renshu.py
# 輸送費の定義
cost = [[10,10,11,17],
        [16,19,12,14],
        [15,12,14,12]]
```

表2に従い、輸送量の変数定義を行います。
輸送量の変数は、それぞれの行と列の座標における変数を"val行_列"として、セルごとに個別に定義します。
iは行位置を示し、1行～3行まであることからそれぞれの値を0,1,2とします。
jは列位置を示し、1列～4列まであることからそれぞれの値を0,1,2,3とします。
輸送量の変数はi,jを用いて以下に定義します。

```7.2.renshu.py
# 輸送量の変数定義
amount = {}
for i in range(3):
    for j in range(4):
        amount[i,j] = solver.IntVar(0, 50, "val%i_%i" % (i,j)) 
```
※注意　"val%i_%i"部分の%iは整数を表しており、(i,j)と表示した部分の"i"とは異なる意味です。

次に制約条件を考えます。
問題文に記載されている制約条件のうち、「各倉庫からの搬出は、供給可能量以下」をまず取り上げます。

倉庫W1における搬出量をamount_w1_sumとします。amount_w1_sumの値は、
`amount_w1_sum = amount[0,0] + amount[0,1] + amount[0,2] + amount[0,3]`
となります。
このamount_w1_sumが供給可能量35を超えないことが制約条件になります。

この条件から、

```7.2.renshu.py
# 制約条件① - 各倉庫からの搬出量は供給可能量以下
amount_w1_sum = 0
for j in range(4):
    amount_w1_sum += amount[0,j]
solver.Add(amount_w1_sum <= 35)
```

倉庫W2,W3においても同様に、

```7.2.renshu.py
amount_w2_sum = 0
for j in range(4):
    amount_w2_sum += amount[1,j]
solver.Add(amount_w2_sum <= 41)
amount_w3_sum = 0
for j in range(4):
    amount_w3_sum += amount[2,j]
solver.Add(amount_w3_sum <= 42)
```

となります。

次に、別の制約条件を考えます。
問題文に記載されている制約条件のうち、「各工場への搬入は、需要量以上」を取り上げます。

工場F1に搬入する搬入量をamount_F1_sumとします。
`amount_F1_sum = amount[0,0] + amount[1,0] + amount[2,0]`
となります。
このamount_F1_sumが需要量28以上となることが制約条件になります。

この条件から、

```7.2.renshu.py
# 制約条件② - 各工場への搬入量は需要量以上
amount_F1_sum = 0
for i in range(3):
    amount_F1_sum += amount[i,0]
solver.Add(amount_F1_sum >= 28)
```

工場F2,F3,F4も同様に、

```7.2.renshu.py
amount_F2_sum = 0
for i in range(3):
    amount_F2_sum += amount[i,1]
solver.Add(amount_F2_sum >= 29)
amount_F3_sum = 0
for i in range(3):
    amount_F3_sum += amount[i,2]
solver.Add(amount_F3_sum >= 31)
amount_F4_sum = 0
for i in range(3):
    amount_F4_sum += amount[i,3]
solver.Add(amount_F4_sum >= 25)
```

本問題の最適化は、「総輸送費が最小となる」計画を求めることです。
最小化問題なので、solver.Minimizeを使用します。
総輸送量は、表1と表2のそれぞれセルごとに輸送量と輸送費をかけたものを全体的に合算することから、最小化の条件は以下となります。

```7.2.renshu.py
# 総輸送費の最小化
solver.Minimize(solver.Sum([amount[i,j] * cost[i][j] for i in range(3) for j in range(4)]))
```

求解を実行します。

```7.2.renshu.py
# 求解実行
status = solver.Solve()
```


結果を確認します。

```7.2.renshu.py
if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    for i in range(3):
        for j in range(4):
            if amount[i,j].SolutionValue() > 0:
                print(i,j,amount[i,j].SolutionValue())
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

最適化計算結果は以下となりました。

>
>Solution:
>ok
>Objective value = 1296.0
>0 0 28.0
>0 1 7.0
>1 2 31.0
>1 3 5.0
>2 1 22.0
>2 3 20.0

以下は実行時間を表示します。
`print("Time = ", solver.WallTime(), " milliseconds")`


プログラム全体です。

```7.2.renshu.py
from __future__ import print_function
from ortools.linear_solver import pywraplp
import numpy as np


# Create the mip solver with the CBC backend.
solver = pywraplp.Solver('simple_mip_program',
    pywraplp.Solver.CBC_MIXED_INTEGER_PROGRAMMING)

# 輸送費の定義
cost = [[10,10,11,17],
        [16,19,12,14],
        [15,12,14,12]]


# 輸送量の変数定義
amount = {}
for i in range(3):
    for j in range(4):
        amount[i,j] = solver.IntVar(0, 50, "val%i_%i" % (i,j)) 


# 制約条件① - 各倉庫からの搬出量は供給可能量以下
amount_w1_sum = 0
for j in range(4):
    amount_w1_sum += amount[0,j]
solver.Add(amount_w1_sum <= 35)
amount_w2_sum = 0
for j in range(4):
    amount_w2_sum += amount[1,j]
solver.Add(amount_w2_sum <= 41)
amount_w3_sum = 0
for j in range(4):
    amount_w3_sum += amount[2,j]
solver.Add(amount_w3_sum <= 42)

# 制約条件② - 各工場への搬入量は需要量以上
amount_F1_sum = 0
for i in range(3):
    amount_F1_sum += amount[i,0]
solver.Add(amount_F1_sum >= 28)
amount_F2_sum = 0
for i in range(3):
    amount_F2_sum += amount[i,1]
solver.Add(amount_F2_sum >= 29)
amount_F3_sum = 0
for i in range(3):
    amount_F3_sum += amount[i,2]
solver.Add(amount_F3_sum >= 31)
amount_F4_sum = 0
for i in range(3):
    amount_F4_sum += amount[i,3]
solver.Add(amount_F4_sum >= 25)


# 総輸送費の最小化
solver.Minimize(solver.Sum([amount[i,j] * cost[i][j] for i in range(3) for j in range(4)]))

# 求解実行
status = solver.Solve()

if status == pywraplp.Solver.OPTIMAL or status == pywraplp.Solver.FEASIBLE:
    print('Solution:')
    print('ok')
    print('Objective value =', solver.Objective().Value())
    for i in range(3):
        for j in range(4):
            if amount[i,j].SolutionValue() > 0:
                print(i,j,amount[i,j].SolutionValue())
    print("Time = ", solver.WallTime(), " milliseconds")
else:
    print('The problem does not have an optimal solution.')
```

#出展

本記事は、数理最適化についての参考テキスト「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」に記載の練習問題を参考にさせて頂いております。

■参考テキスト
「Pythonによる問題解決シリーズ　データ分析ライブラリーを用いた最適化モデルの作り方」
　斉藤努[著]
　近代科学社[出版]

![001.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/399135/df8e6661-9dd6-559b-5ae8-0bf79f73d1a5.jpeg)

# ご意見など

ご意見、間違い訂正などございましたらお寄せ下さい。