Python + PuLP でたくさん輸送問題
あれでもいいのですが、あのコードでは一回解いたきりで終わってしまいます。
ただの趣味ならいいけれど、最適化を外部から呼び出せるものとして実装する必要が(普通は)あるでしょう。そして、ランダムでインスタンス(問題例)が作れるタイプの問題であれば、ガンガンインスタンスを作ってガンガン解くようなテスターを作っておくと色々と便利です。定式化や解法の違いによる計算時間の比較や問題サイズに対する計算量の増加も見やすくなります。
私は、カプセル化のようなことをして、モデルオブジェクトと変数オブジェクトとその他必要なデータとか出力用の関数とかを詰め込んだラッパ (wrapperです、念のため) のクラスをよく作ります。このとき,同じファイルのメイン関数で1個だけインスタンス(問題例)を解くような簡易テスト付き.pyファイルに仕立ててあげるとデバッグも楽になります。ランダムなインスタンスを作って解くテスターは別のpyファイルとかJupyter Notebookに書き、ラッパをimportして使うとすっきりします。
プログラミングの腕も開発経験もない私のやり方なので、他のやり方もあると思います。ひとつの参考程度になさってください。
輸送問題のラッパ
実行すると簡単な問題例を解いて結果を表示します。
""" transportation.py Local Optima 2020. """ import itertools as itrt import pulp EPS = 1e-4 def _main(): capacity = dict(enumerate([6,10])) demand = dict(enumerate([4, 6, 3, 3])) transcost = {(0, 0): 3, (0, 1): 7, (0, 2): 11, (0, 3): 8, (1, 0): 6, (1, 1): 7, (1, 2): 8, (1, 3): 9} tp = Transportation(capacity, demand, transcost) tp.solve() tp.print_solution() class Transportation(): def __init__(self, capacity, demand, transcost): """ Transportation Problem Interface Parameters----- capacity: dict of factory capacity demand: dict of demand qty from customers transcost: dict of transportation cost from i to j whose keys are tuple (i,j) """ self.capacity = capacity self.demand = demand self.transcost = transcost def solve(self): """ Solve transportation problem """ model, x = self._formulate() model.solve() print("Transportation Problem Solved --", pulp.LpStatus[model.status]) def _formulate(self): model = pulp.LpProblem("Transportation", pulp.LpMinimize) # 変数オブジェクト x = {(i,j): pulp.LpVariable("x{}-{}".format(i,j), lowBound=0, upBound=self.demand[j]) for i,j in self.transcost} # 目的関数 model += pulp.lpSum([self.transcost[i,j] * x[i,j] for i,j in self.transcost]), "Objective" # 制約条件 # 工場の出荷上限 for i, Ci in self.capacity.items(): model += pulp.lpSum([x[i,j] for j in self.demand.keys()]) <= Ci, "Capacity{}".format(i) # 顧客の需要 for j, dj in self.demand.items(): model += pulp.lpSum([x[i,j] for i in self.capacity.keys()]) == dj, "demand{}".format(j) self.model = model self.x = x return model, x def print_solution(self): """ Print Solution if possible """ if self.model.status != 1: print("Not Ready to print.") return print("z =", pulp.value(self.model.objective)) for var in self.model.variables(): if var.varValue > EPS: print(var, var.varValue) if __name__ == "__main__": _main()
テスター
問題例をたくさん生成して解きまくるテスターは,たとえば以下のようになります。
さきほどのtransportation.pyをfrom transportation import *
して使うのがミソです。
import math import random import itertools as itrt from transportation import * def distance(pos1, pos2): return sum([(x1-x2)**2 for x1,x2 in zip(pos1, pos2)])**0.5 def random_transportation(factory, store, mapsize=10): """ Returns Transportation problem wrapper Parameters----- factory: number of factories store: number of stores mapsize (optional): """ capacity = dict(enumerate([random.randint(3,15) for _ in range(factory)])) demand = dict(enumerate([random.randint(1,10) for _ in range(store)])) posi = {n:[random.uniform(0,mapsize), random.uniform(0,10)] for n in range(factory)} posj = {n:[random.uniform(0,mapsize), random.uniform(0,10)] for n in range(store)} dist = {} for i,j in itrt.product(posi, posj): dist[i,j] = math.ceil(distance(posi[i], posj[j])) return Transportation(capacity, demand, dist) for _ in range(30): tp = random_transportation(4,6, mapsize=20) tp.solve()
実行結果です。
Transportation Problem Solved -- Optimal Transportation Problem Solved -- Infeasible Transportation Problem Solved -- Optimal Transportation Problem Solved -- Optimal Transportation Problem Solved -- Optimal (以下省略)
capacity
とdemand
をランダムに作っているので,工場の総供給能力が需要の総合計に負けてしまう場合は実行不可能(Infeasible)になります。
補遺
輸送問題は,実は面白い性質をもった線形計画問題です。
各工場の供給能力と,各店舗の需要量を整数値にしておくと,なんと連続変数なのに最適解が整数になるのです。気になる人は上のテスターをいじくって試してみてください。
これはたまたまではないことが知られています。
輸送問題を の形に書いたときの行列 のことを制約行列 といいますが,これが完全単模 (totally-unimodular) であるとき,連続LPを単体法で解いたときの最適解が整数になるそうです。