PuLPで最適化シリーズ第2回は、輸送問題を解いてみます。工場から店舗への供給を最適化する、ロジスティクスの基本問題です。

モデリングの復習でもありますが，リストや辞書を使ったPythonらしい書き方でモデリングをしていきます。

輸送問題の定義

ｍ個の工場 $0,\cdots,i,\cdots,m-1$ から，n個の店舗 $0,\cdots,j,\cdots,n-1$ に物資を輸送します。工場 $i$ の供給上限 $C_i$ 、店舗jの需要量 $d_j$ を満たしつつ、輸送にかかる総費用が最小となるような、工場iから店舗jへの出荷量 $x _ {ij}$ を求めよ。工場iから店舗jへの単位輸送費用は $c _ {ij}$ である。

輸送問題では簡単のため，総輸送費用を

$\displaystyle{\sum _ {0 \le i \lt m}\sum _ {0 \le j \lt n} c _ {ij} x _ {ij}}$

とします。

（Pythonのリストを使ってこのあとの処理をする都合上、添字をゼロスタートで考えてます。）

このとき輸送問題は次の最小化の線形計画問題になります。

$\begin{aligned} \text{minimize}\quad & \sum _ {0 \le i \lt m}\sum _ {0 \le j \lt n} c _ {ij} x _ {ij} \\ \text{subject to}\quad & \sum _ {0 \le j \lt n} x _ {ij} \le C_i & (0 \le i \lt n) \\ & \sum _ {0 \le i \lt m} x _ {ij} = d_j & (0 \le j \lt m) \\ & \forall x _ {ij } \ge 0 \end{aligned}$

目的関数は総費用の最小化。

制約条件の1つ目は工場の出荷量合計が供給上限を超えない制約，2つ目は店舗の入荷量合計が需要と一致する制約です。

例題

Pythonのリストで例題のデータを作ります。 $m=2, n=4$ の例を見てみましょう。

その前に，Python + PuLPで最適化するときの段取りも思い出してくださいね。

インポートする (from pulp import *)
モデルオブジェクトを宣言する
変数オブジェクトを宣言する
目的関数をセットする
制約条件をセットする
最適化する
結果を見る

from pulp import *

# 工場のリスト（中身は供給上限）
capacity = [6, 10]

# 店舗のリスト（中身は需要量）
demand = [4, 6, 3, 3]

# 単位費用の行列 c_ij
cost = {(0, 0): 3,
        (0, 1): 7,
        (0, 2): 11,
        (0, 3): 8,
        (1, 0): 6,
        (1, 1): 7,
        (1, 2): 8,
        (1, 3): 9  }

工場と店舗のリストはこういう風にリストで与えるといいです。

単位費用の行列 $(c _ {ij} )$ は，二重リストにすることもできますが，私のオススメはタプルをキーとする辞書(dict)にする方法です。添字のたくさんつくものは辞書にせよ。もし今後 $c _ {ij} ^k$ みたいなものがモデルに現れても，cost[i,j,k]で参照できるようなdictにするとよいです。

なぜか。Zen of Pythonにも書いてあるとおり，"Flat is better than nested." だからです。でも，どうしても多重リストが好きな人は多重リストにしてください。

次は，モデルオブジェクトと変数オブジェクトの宣言です。変数も2添字の行列の形になりますので，for文の中で回すことになります。

model = LpProblem("Transportation", LpMinimize)

# 変数オブジェクト
x = {}
for i,j in cost:
    x[i,j] = LpVariable("x{}-{}".format(i,j), lowBound=0, upBound=demand[j])

工場iから店舗jへの出荷量 $x _ {ij}$ は， $+\infty$ にデカくなることはなく，高々，店舗jの需要量までなのは間違いないので，upBound=demand[j]を入れてあげますが，敢えて入れなくても問題ありません。

リスト内包表記に強い人は，変数のところを

x = {(i,j): LpVariable("x{}-{}".format(i,j), lowBound=0, upBound=demand[j]) for i,j in cost}

と1行にしてもかまいません。

目的関数と制約条件の設定です。和の形になっているものは，ベタ書きしないでも，LpSumという便利な関数があります。 lpSum([リスト内包表記])の形で使ってあげます。

# 目的関数
model += lpSum([cost[i,j] * x[i,j] for i,j in cost]), "Objective"

# 制約条件
# 工場の出荷上限
for i, Ci in enumerate(capacity):
    model += lpSum([x[i,j] for j in range(len(demand))]) <= Ci, "Capacity{}".format(i)

# 顧客の需要
for j, dj in enumerate(demand):
    model += lpSum([x[i,j] for i in range(len(capacity))]) == dj, "demand{}".format(j)

あとは最適化して結果を見てみましょう。

# 求解
model.solve()

# 結果の確認
print("Optimal Value =", value(model.objective))
for var in x.values():
    # 誤差以上の値を持っている変数だけprint
    if var.varValue > 1e-4:
        print(var, var.varValue)

ソルバーの挙動として，浮動小数点に誤差が乗って解となる場合があります。そのため，値1e-4 $=10^{-4}$ 以上の変数だけ表示してあげています。

Optimal Value = 103.0
x0_0 4.0
x0_3 2.0
x1_1 6.0
x1_2 3.0
x1_3 1.0

こんな結果でした。店舗3には２つの工場から出荷することになっていますね。

データの生成

距離行列は，適当に作るとウソくさいデータセットになるので，

10x10のマップにランダムに工場と店舗を配置し
2地点間のユークリッド距離を切り上げ

た距離行列にしてみました。以下のコードで作れます。

import math
import random
import itertools as itrt
from pprint import pprint

def distance(pos1, pos2):
    return sum([(x1-x2)**2 for x1,x2 in zip(pos1, pos2)])**0.5

posi = {n:[random.uniform(0,10), random.uniform(0,10)] for n in range(2)}
posj = {n:[random.uniform(0,10), random.uniform(0,10)] for n in range(4)}

dist = {}
for i,j in itrt.product(posi, posj):
    dist[i,j] = math.ceil(distance(posi[i], posj[j]))
pprint(dist)

itertoolsは私はas itrtってインポートしてるんですが，みなさんどうしてますか？よかったらコメント欄で教えて下さい。

Local Optima

数理最適化を趣味的に楽しむ。

Python + PuLPで輸送問題

輸送問題の定義

例題

データの生成