Local Optima

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

Python + PuLP でたくさん輸送問題

第2回で輸送問題をPuLPで解きました。

local-optima.hateblo.jp

あれでもいいのですが、あのコードでは一回解いたきりで終わってしまいます。

ただの趣味ならいいけれど、最適化を外部から呼び出せるものとして実装する必要が(普通は)あるでしょう。そして、ランダムでインスタンス(問題例)が作れるタイプの問題であれば、ガンガンインスタンスを作ってガンガン解くようなテスターを作っておくと色々と便利です。定式化や解法の違いによる計算時間の比較や問題サイズに対する計算量の増加も見やすくなります。

 

私は、カプセル化のようなことをして、モデルオブジェクトと変数オブジェクトとその他必要なデータとか出力用の関数とかを詰め込んだラッパ (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
(以下省略)

capacitydemandをランダムに作っているので,工場の総供給能力が需要の総合計に負けてしまう場合は実行不可能(Infeasible)になります。

補遺

輸送問題は,実は面白い性質をもった線形計画問題です。

各工場の供給能力と,各店舗の需要量を整数値にしておくと,なんと連続変数なのに最適解が整数になるのです。気になる人は上のテスターをいじくって試してみてください。

これはたまたまではないことが知られています。

輸送問題を  \text{min} \{ cx | Ax \ge b \} の形に書いたときの行列  A のことを制約行列 といいますが,これが完全単模 (totally-unimodular) であるとき,連続LPを単体法で解いたときの最適解が整数になるそうです。