import random
import math
from gurobipy import Model, quicksum, GRB, tuplelist, multidict
# from mypulp import Model, quicksum, GRB, tuplelist, multidict
import networkx as nx
import matplotlib.pyplot as plt
import plotly
多品種流問題
最小費用流問題の拡張として， 複数の異なる「もの」のフローを扱う多品種流問題（multi-commodity flow problem）を考えよう．
ネットワーク上を流す必要がある異なる「もの」を品種（commodity）とよぶ． 品種は，始点と終点をもち，決められた量を始点から終点まで運ぶ必要があるものとする． さて，これらの品種に対して，別々に容量制約が与えられている場合には， 品種ごとに前節の最小費用流問題を解けば済む． ここでは，自然な拡張として，品種の流量の和が， 枝の容量以下であるという条件を課すものとする．
多品種流問題
$n$ 個の点から構成される点集合 $V$ および $m$ 本の枝から構成される枝集合 $E$, 品種集合 $K$， $V$ および $E$ から成る有向グラフ $G=(V,E)$， 品種ごとに枝上に定義されるフロー1単位あたりの費用関数 $c: E \times K \rightarrow \mathbf{R}$， 枝上に定義される非負の容量関数 $u: E \rightarrow \mathbf{R}_+ \cup \{ \infty \}$， 品種 $k(\in K)$ の始点 $s_k$，終点 $t_k$，品種の需要量関数 $b: K \rightarrow \mathbf{R}$ が与えられたとき， 「実行可能フロー」で，費用の合計が最小になるものを求めよ．
多品種流問題における実行可能フローとは，フロー整合条件，容量制約，ならびに非負制約を満たすフローである． 正確に言うと，多品種の実行可能フローとは， 実数値関数 $x: E \times K \rightarrow \mathbf{R}$ で，以下の性質を満たすものを指す．
フロー整合条件: $$ \sum_{j: ji \in E} x_{ji}^k - \sum_{j: ij \in E} x_{ij}^k = \left\{ \begin{array}{ll} -b_k & i=s_k \\ 0 & \forall i \in V \setminus \{s_k,t_k\} \\ b_k & i=t_k \end{array} \right. ,\forall k \in K $$
容量制約: $$ \sum_{k \in K} x_{e}^k \leq u_{e} \ \ \ \forall e \in E $$
非負制約: $$ x_{e}^k \geq 0 \ \ \ \forall e \in E, k \in K $$
上の定義では，品種は分割して運んでも良いものと仮定した． このことを強調するために，上の問題は小数多品種流問題（fractional multi-commodity flow problem）と よばれることもある．小数多品種流問題は，線形最適化問題の特殊形であるので，多項式時間で解くことができる．
一方，分割を許さない多品種フロー問題は，整数多品種流問題（integer multi-commodity flow problem） とよばれ，品種数が $2$ で容量がすべて $1$ の場合でも$NP$-困難になることが知られている．
多品種流問題は，無向グラフ $G=(V,E)$ でも自然に定義できる． 無向グラフの場合には，枝 $e=ij$ 上の容量制約は以下のように変更される．
$$ \sum_{k \in K} \left( x_{ij}^k + x_{ji}^k \right) \leq u_{e} \ \ \ \forall e \in E $$これは，$i$ から $j$ に流れるフローと $j$ から $i$ に流れるフローの和が，枝 $ij$ の容量を超えないことを表す．
上の制約の下で，費用の合計 $ \sum_{k \in K} \sum_{e \in E} c_e^k x_{e}^k $ を最小化する．
m, n = 3, 3
cost_lb, cost_ub = 10, 10
cap_lb, cap_ub = 150, 150
demand_lb, demand_ub = 10, 30
G = nx.grid_2d_graph(m, n)
D = G.to_directed()
for (i, j) in D.edges():
D[i][j]["cost"] = random.randint(cost_lb, cost_ub)
D[i][j]["capacity"] = random.randint(cap_lb, cap_ub)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
b = {}
K = []
for i in D.nodes():
for j in D.nodes():
if i != j:
K.append((i, j))
b[i, j] = random.randint(demand_lb, demand_ub)
V = set(D.nodes())
E = set(D.edges())
model = Model("multi-commodity flow")
x = {}
for k in K:
for (i, j) in E:
x[i, j, k] = model.addVar(vtype="C", name=f"x(i,j,k)")
model.update()
for i in V:
for k in K:
if i == k[0]:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== -b[k]
)
elif i == k[1]:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== b[k]
)
else:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== 0
)
# Capacity constraints
for (i, j) in E:
model.addConstr(
quicksum(x[i, j, k] for k in K) <= D[i][j]["capacity"], f"Capacity({i},{j})"
)
# Objective:
model.setObjective(
quicksum(D[i][j]["cost"] * x[i, j, k] for (i, j, k) in x), GRB.MINIMIZE
)
model.optimize()
flow = {}
for (i, j) in E:
flow[i, j] = sum(x[i, j, k].X for k in K)
plt.figure()
nx.draw(D, pos=pos, with_labels=True, node_size=1000, node_color="yellow")
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{flow[i,j]} \n{flow[j,i]}"
nx.draw_networkx_edge_labels(D, pos, edge_labels=edge_labels)
plt.show()
多品種輸送問題
今度は，複数の製品（これを品種と呼ぶ）を運ぶための輸送問題の拡張を考える．
いま，複数の製品を工場から顧客へ輸送することを考える． ただし，製品によって重量が異なるため輸送費用が変わってくるものと仮定する． この問題は，多品種輸送問題(multi-commodity transportation problem)と呼ばれる．
定式化に必要な記号は，輸送問題とほとんど同じであるので，相違点のみ定義しておく． 製品（品種）の集合を $K$ とする． 工場 $j$ から顧客 $i$ に製品 $k$ が輸送される量を表す連続変数 $x_{ijk}$ を導入する．
顧客 $i$ における製品 $k$ の需要量を $d_{ik}$， 工場 $j$ の生産量上限（容量）を $M_j$， 顧客 $i$ と工場 $j$ 間に製品 $k$ の $1$ 単位の需要が移動するときにかかる 輸送費用を $c_{ijk}$ とする．
上の記号を用いると，多品種輸送問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \displaystyle\sum_{i \in I} \displaystyle\sum_{j \in J} \displaystyle\sum_{k \in K} c_{ijk} x_{ijk} & \\ s.t. & \displaystyle\sum_{j \in J} x_{ijk} =d_{ik} & \forall i \in I, k \in K \\ & \displaystyle\sum_{i \in I} \displaystyle\sum_{k \in K} x_{ijk} \leq M_j & \forall j \in J \\ & x_{ijk} \geq 0 & \forall i \in I, j \in J, k \in K \end{array} $$最初の制約は，各製品ごとに需要が満たされることを表し， 2番目の制約は，工場で生産されるすべての製品の合計量が，工場の容量を超えないことを表す．
def mctransp(I, J, K, c, d, M):
"""mctransp -- model for solving the Multi-commodity Transportation Problem
Parameters:
- I: set of customers
- J: set of facilities
- K: set of commodities
- c[i,j,k]: unit transportation cost on arc (i,j) for commodity k
- d[i][k]: demand for commodity k at node i
- M[j]: capacity
Returns a model, ready to be solved.
"""
model = Model("multi-commodity transportation")
# Create variables
x = {}
for (i, j, k) in c:
x[i, j, k] = model.addVar(vtype="C", name=f"x(i,j,k)")
model.update()
arcs = tuplelist([(i, j, k) for (i, j, k) in x])
# Demand constraints
for i in I:
for k in K:
model.addConstr(
quicksum(x[i, j, k] for (i, j, k) in arcs.select(i, "*", k)) == d[i, k],
f"Demand(i,k)",
)
# Capacity constraints
for j in J:
model.addConstr(
quicksum(x[i, j, k] for (i, j, k) in arcs.select("*", j, "*")) <= M[j],
f"Capacity(j)",
)
# Objective:
model.setObjective(
quicksum(c[i, j, k] * x[i, j, k] for (i, j, k) in x), GRB.MINIMIZE
)
model.update()
model.__data = x
return model
def make_inst1():
d = {
(1, 1): 80,
(1, 2): 85,
(1, 3): 300,
(1, 4): 6, # {customer: {commodity:demand <float>}}
(2, 1): 270,
(2, 2): 160,
(2, 3): 400,
(2, 4): 7,
(3, 1): 250,
(3, 2): 130,
(3, 3): 350,
(3, 4): 4,
(4, 1): 160,
(4, 2): 60,
(4, 3): 200,
(4, 4): 3,
(5, 1): 180,
(5, 2): 40,
(5, 3): 150,
(5, 4): 5,
}
I = set([i for (i, k) in d])
K = set([k for (i, k) in d])
J, M = multidict({1: 3000, 2: 3000, 3: 3000}) # capacity
produce = {
1: [2, 4],
2: [1, 2, 3],
3: [2, 3, 4],
} # products that can be produced in each facility
weight = {1: 5, 2: 2, 3: 3, 4: 4} # {commodity: weight<float>}
cost = {
(1, 1): 4,
(1, 2): 6,
(1, 3): 9, # {(customer,factory): cost<float>}
(2, 1): 5,
(2, 2): 4,
(2, 3): 7,
(3, 1): 6,
(3, 2): 3,
(3, 3): 4,
(4, 1): 8,
(4, 2): 5,
(4, 3): 3,
(5, 1): 10,
(5, 2): 8,
(5, 3): 4,
}
c = {}
for i in I:
for j in J:
for k in produce[j]:
c[i, j, k] = cost[i, j] * weight[k]
return I, J, K, c, d, M
I, J, K, c, d, M = make_inst1()
model = mctransp(I, J, K, c, d, M)
model.optimize()
print("Optimal value:", model.ObjVal)
EPS = 1.0e-6
x = model.__data
for (i, j, k) in x:
if x[i, j, k].X > EPS:
print(
"sending %10s units of %3s from plant %3s to customer %3s"
% (x[i, j, k].X, k, j, i)
)
多品種ネットワーク設計問題
多品種流問題に枝上の固定費用 $F: E \rightarrow \mathbf{R}_+$ をつけた問題を，多品種ネットワーク設計問題(multi-commodity network design problem)とよぶ． この問題は，$NP$-困難であり，しかも以下の通常の定式化だと小規模の問題例しか解けない．実務的には，パス型の定式化や列生成法を用いることが推奨される．
枝を使用するか否かを表す $0$-$1$変数を用いる．
目的関数： $$ \sum_{k \in K} \sum_{e \in E} c_e^k x_{e}^k + \sum_{e \in E} F_{ij} y_{ij} $$
フロー整合条件： $$ \sum_{j: ji \in E} x_{ji}^k - \sum_{j: ij \in E} x_{ij}^k = \left\{ \begin{array}{ll} -b_k & i=s_k \\ 0 & \forall i \in V \setminus \{s_k,t_k\} \\ b_k & i=t_k \end{array} \right. ,\forall k \in K $$
容量制約: $$ \sum_{k \in K} x_{ij}^k \leq u_{ij} y_{ij} \ \ \ \forall (i,j) \in E $$
非負制約: $$ x_{e}^k \geq 0 \ \ \ \forall e \in E, k \in K $$
整数制約： $$ y_{ij} \in \{ 0,1 \} \ \ \ \forall (i,j) \in E $$
m, n = 4, 4
cost_lb, cost_ub = 10, 10
fixed_lb, fixed_ub = 1000, 1000
cap_lb, cap_ub = 1500, 1500
demand_lb, demand_ub = 10, 30
G = nx.grid_2d_graph(m, n)
D = G.to_directed()
for (i, j) in D.edges():
D[i][j]["cost"] = random.randint(cost_lb, cost_ub)
D[i][j]["fixed"] = random.randint(fixed_lb, fixed_ub)
D[i][j]["capacity"] = random.randint(cap_lb, cap_ub)
pos = {(i, j): (i, j) for (i, j) in G.nodes()}
b = {}
K = []
for i in D.nodes():
for j in D.nodes():
if i != j:
K.append((i, j))
b[i, j] = random.randint(demand_lb, demand_ub)
V = set(D.nodes())
E = set(D.edges())
model = Model("multi-commodity design")
x, y = {}, {}
for (i, j) in E:
y[i, j] = model.addVar(vtype="B", name=f"y(i,j)")
for k in K:
x[i, j, k] = model.addVar(vtype="C", name=f"x(i,j,k)")
model.update()
for i in V:
for k in K:
if i == k[0]:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== -b[k]
)
elif i == k[1]:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== b[k]
)
else:
model.addConstr(
quicksum(x[j, i, k] for j in V - {i} if (j, i) in E)
- quicksum(x[i, j, k] for j in V - {i} if (i, j) in E)
== 0
)
# Capacity constraints
for (i, j) in E:
model.addConstr(
quicksum(x[i, j, k] for k in K) <= D[i][j]["capacity"] * y[i, j],
f"Capacity({i},{j})",
)
# Objective:
model.setObjective(
quicksum(D[i][j]["cost"] * x[i, j, k] for (i, j, k) in x)
+ quicksum(D[i][j]["fixed"] * y[i, j] for (i, j) in y),
GRB.MINIMIZE,
)
model.optimize()
flow = {}
for (i, j) in E:
flow[i, j] = int(sum(x[i, j, k].X for k in K))
plt.figure()
nx.draw(D, pos=pos, with_labels=True, node_size=1000, node_color="yellow")
edge_labels = {}
for (i, j) in G.edges():
edge_labels[i, j] = f"{flow[i,j]} \n{flow[j,i]}"
nx.draw_networkx_edge_labels(D, pos, edge_labels=edge_labels)
path = []
for (i, j, k) in x:
if k == ((0, 0), (n - 1, n - 1)): # (0,0)から(3,3)へのパスを表示
if x[i, j, k].X > 0.001:
# print(i,j,k,x[i,j,k].X)
path.append((i, j))
nx.draw_networkx_edges(G, pos, edgelist=path, width=10, edge_color="orange")
plt.show()
サービス・ネットワーク設計問題
ロジスティクス・ネットワークが供給地点から需要地点への，一方向的な「もの」の流れを扱うのに対して， サービス・ネットワークでは，発地と着地の間に輸送される，多対多の「もの」の流れを扱う． ロジスティクス・ネットワークにおいては，顧客需要の不確実性やサービス・レベルの要求などの条件を満たすために， ネットワーク内の地点で在庫を適切に管理することが重要になる． それに対して，サービス・ネットワークにおいては，基本的には途中で在庫することなく， 発地から着地へ「もの」が流れていく． （ただし，輸送のタイミングをとるために，一時保管することは許される．）
サービス・ネットワーク設計問題（consolidation transportation network design problem）は， しばしば混載ネットワーク設計問題ともよばれる． しかし，我が国における実務家の間では，「混載」という用語が「異なる種類の製品の積み合わせ」という意味で用いられ，本モデルでは，同じ種類の製品の積み合わせが主な応用であるため， 混乱を避ける意味でサービス・ネットワーク設計問題とよぶことにする． ここで考えるモデルでは， 主に（郵便物や宅配便のように）発地・着地が異なる同じ種類の製品の積み合わせを対象とする 場合が多いが，異なる種類の製品の積み合わせ（いわゆる業界用語での混載）への適用も可能である．
配送計画は，比較的短距離の輸送を計画するためのモデルであるが， サービス・ネットワーク設計では比較的長距離の輸送を対象とする． そのため，途中での積み替えを考慮する必要が出てくる． これが，ここで考えるサービス・ネットワーク設計が， 配送計画と異なったアプローチを必要とする理由である． 配送計画においては，積み替えを考慮する必要がなかったので，運搬車の移動経路を求めれば， 荷（品目）の移動は自動的に定まった．一方，サービス・ネットワーク設計においては， 積み替えを考慮する必要があるので，運搬車の流れのみならず，荷の流れが意思決定項目となり， これによって問題の難しさが増大する．
ここで考えるモデルは，実務から生まれたものであり，以下の仮定に基づく．
荷（load）とは，地点間を移動させる「モノ」の総称である．ネットワーク理論では品種(commodity) とよばれるが，ここでは現実問題を想定していることを強調するために荷とよぶことにする． 荷は，その移動によって利益を生む資源の総称である． 荷が最初に出発する（最後に到着する）地点を荷の発地（着地）とよび，あらかじめ決められている． ネットワーク理論では，発地（source）は始点もしくは発生点，着地（sink, terminal）は終点もしくは集中点とよばれるが，以下では，発地，着地とよぶものとする．
荷は途中で分岐してはならない．言い換えれば，荷の発地から着地までの移動経路は，1本のパスでなければならない．
宅配便や郵便事業などへの応用では，さらに荷の移動に入木条件とよばれる制約が付加される． 入木条件とは，着地が同じ荷のパスが合流したら，その後のパスは同じ経路にならなければならないことを規定する． これは，集荷した荷物を積み替え（ソーティング）する際に，現場でのオペレーションの簡便性のため， 着地の情報だけを利用するためである． 実際に，東京行きの荷物という荷の集まりは，様々な発地からの荷が集約されたものであり，これを発地別に分けて 異なる方面行きの運搬車に積み替えることは （たとえそれによって費用が減少する可能性があっても）現実には行われないのである．
運搬車（トラック）の積載容量は既知であり，積載される荷量の合計は積載容量を超えない．
地点間の運搬車の移動費用は既知である．
中継地点での積替え費用は荷量によって定まり，既知である．
運搬車は出発地点に戻ってくる必要はない．（この仮定を緩めるのは比較的容易である．）
上のサービス・ネットワーク設計モデルに対して，数理最適化（ならびに制約最適化）ソルバーを利用した専用解法が構築できる． 専用解法を組み込んだサービス・ネットワーク最適化システムとしてSENDO (https://www.logopt.com/sendo) が開発されている．