blog

最小費用流問題

Published:

By nob

Category: Posts

Tags: 数理最適化 Python-MIP Python

燃料費高騰のため輸送費を削減したい

工場で製造した製品を需要地に近い倉庫に車両で運搬する。燃料費高騰のため輸送費を削減したい。

  • 工場の出荷可能量 
supply = [39, 40, 41]
  • 倉庫の需要量 
demand = [38, 20, 22, 36]
  • 輸送費の単価（1商品あたり） 
cost = [
    [19, 12, 13, 18],
    [14, 12, 18, 12],
    [15, 16, 15, 10],
]

モデルの作成 

from mip import Model, OptimizationStatus, minimize, xsum

model = Model(solver_name="cbc")
model.verbose = 0

product = model.add_var_tensor((len(supply), len(demand)), "product")
model.objective = minimize(xsum((cost * product).flat))

for f in range(len(supply)):
    model += xsum(product[f]) <= supply[f]

for w in range(len(demand)):
    model += xsum(product[:, w]) == demand[w]

model.optimize()

if model.status == OptimizationStatus.OPTIMAL:
    print(model.objective_value)
    print(product.astype(float, subok=False))
else:
    print(model.status)

1419.0
[[ 0. 17. 22.  0.]
 [37.  3.  0.  0.]
 [ 1.  0.  0. 36.]]

グラフの表示 

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np


def draw_graph(supply, demand, solution):

    fig, ax = plt.subplots(1, 1)

    g = nx.DiGraph()

    pos = {}

    for i, s in enumerate(supply):
        node_name = "F" + str(i + 1)
        g.add_node(node_name, amount=s)
        pos[node_name] = (0, i + 0.5)

    for i, d in enumerate(demand):
        node_name = "W" + str(i + 1)
        g.add_node(node_name, amount=d)
        pos[node_name] = (1, i)

    for f, w in zip(*np.where(solution > 0)):
        g.add_edge("F" + str(f + 1), "W" + str(w + 1), product=solution[f, w])

    node_size = list(
        map(
            lambda a: 50 * a,
            nx.get_node_attributes(g, "amount").values(),
        )
    )

    nx.draw_networkx_nodes(
        g,
        pos,
        node_size=node_size,
        node_color="white",
        edgecolors="black",
        linewidths=3,
        ax=ax,
    )

    edge_width = list(
        map(
            lambda p: 0.5 * p,
            nx.get_edge_attributes(g, "product").values(),
        )
    )

    nx.draw_networkx_edges(
        g,
        pos,
        width=edge_width,
        node_size=1000,
        edge_color="black",
        arrowsize=40,
        ax=ax,
    )

    nx.draw_networkx_labels(g, pos, font_size=15)

    edge_label = nx.get_edge_attributes(g, "product")
    nx.draw_networkx_edge_labels(
        g, pos, edge_label, label_pos=0.2, font_size=12,
        ax=ax,
    )

    plt.axis("off")
    plt.tight_layout()
    plt.show()

draw_graph(supply, demand, product.astype(float, subok=False))

fig-1

積載率が低い！

車両の容積(capacity)を確認したところ 

capacity = 5

だった。

F3 -> W1 は1個、F3 -> W4 は36個運搬するので、必要便数は

  • F3 -> W1: math.ceil(1 / 5) = 1

  • F3 -> W4: math.ceil(36 / 5) = 8

出荷先倉庫によって車両の積載率が大きく異なる。積載率が80%以上になるようにしたい。

1便あたりの単価は積載量にはあまり影響を受けないと考えてモデルを再作成する。

  • 工場の出荷可能量 
supply = [39, 40, 45]
  • 倉庫の需要量 
demand = [38, 20, 22, 36]
  • 輸送費の単価（1便あたり） 
cost_per_transport = np.ceil(np.array(cost) * capacity * 0.9)
cost_per_transport

array([[86., 54., 59., 81.],
       [63., 54., 81., 54.],
       [68., 72., 68., 45.]])

モデルの再作成 

model = Model(solver_name="cbc")
model.verbose = 0

product = model.add_var_tensor((len(supply), len(demand)), "product")
transport = model.add_var_tensor(product.shape, "transport", var_type="I")

model.objective = minimize(xsum((cost_per_transport * transport).flat))

for f in range(len(supply)):
    for w in range(len(demand)):
        model += transport[f, w] >= product[f, w] / capacity
        model += product[f, w] >= 0.8 * transport[f, w] * capacity

for f in range(len(supply)):
    model += xsum(product[f]) <= supply[f]

for w in range(len(demand)):
    model += xsum(product[:, w]) == demand[w]

model.optimize()

if model.status == OptimizationStatus.OPTIMAL:
    print(model.objective_value)
    print(product.astype(float, subok=False))
    print(transport.astype(float, subok=False))
else:
    print(model.status)

1380.0
[[ 0. 15. 22.  0.]
 [34.  5.  0.  0.]
 [ 4.  0.  0. 36.]]
[[0. 3. 5. 0.]
 [7. 1. 0. 0.]
 [1. 0. 0. 8.]]

積載率の計算 

p = product.astype(float, subok=False)
t = transport.astype(float, subok=False)
t = np.where(t > 0, t, 1)
p / t / capacity

array([[0.        , 1.        , 0.88      , 0.        ],
       [0.97142857, 1.        , 0.        , 0.        ],
       [0.8       , 0.        , 0.        , 0.9       ]])

グラフの表示 

draw_graph(supply, demand, product.astype(float, subok=False))

fig-2

データの出典