k-メディアン問題に対する定式化
%matplotlib inline
import random
import math
from gurobipy import Model, quicksum, GRB, multidict
# from mypulp import Model, quicksum, GRB, multidict
import networkx as nx
import matplotlib.pyplot as plt
from amplpy import AMPL, Environment, DataFrame
AMPL_PATH = r"/Users/mikiokubo/Dropbox/My Mac (mikionoiMac-Pro.local)/Documents/ampl/ampl.macos64"
メディアン問題
メディアン問題は，顧客から最も近い施設への距離の「合計」を最小にするように グラフ内の点上または空間内の任意の点から施設を選択する施設配置問題の総称である．
メディアン問題においては，選択される施設の数があらかじめ決められていることが多く， その場合には選択する施設数 $k$ を頭につけて $k$-メディアン問題($k$-median problem)とよばれる． 施設数を表す記号としては基本的にはどんな文字でも良いが， 慣用では $p$ または $k$ を用いることが多いようである． 以下では $k$ を用いることにする．
顧客 $i$ から施設 $j$ への距離を $c_{ij}$ とし，以下の変数を導入する．
$$ x_{ij}= \left\{ \begin{array}{ll} 1 & 顧客 i の需要が施設 j によって満たされるとき \\ 0 & それ以外のとき \end{array} \right. $$$$ y_j = \left\{ \begin{array}{ll} 1 & 施設 j を開設するとき \\ 0 & それ以外のとき \end{array} \right. $$顧客数を $n$ とし，顧客の集合を $I$，施設の配置可能な点の集合を $J$ とする． 通常の $k$-メディアン問題では，施設の候補地点は顧客上と仮定するため，$I=J=\{1,2,\ldots,n\}$ となる．
上の記号および変数を用いると，$k$-メディアン問題は以下の整数最適化問題 として定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \sum_{j \in J} x_{ij} =1 & \forall i \in I \\ & \sum_{j \in J} y_{j} = k & \\ & x_{ij} \leq y_j & \forall i \in I, j \in J \\ & x_{ij} \in \{ 0,1 \} & \forall i \in I, j \in J \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array} $$def kmedian(I, J, c, k):
"""kmedian -- minimize total cost of servicing
customers from k facilities
Parameters:
- I: set of customers
- J: set of potential facilities
- c[i,j]: cost of servicing customer i from facility j
- k: number of facilities to be used
Returns a model, ready to be solved.
"""
model = Model("k-median")
x, y = {}, {}
for j in J:
y[j] = model.addVar(vtype="B", name="y(%s)" % j)
for i in I:
x[i, j] = model.addVar(vtype="B", name="x(%s,%s)" % (i, j))
model.update()
for i in I:
model.addConstr(quicksum(x[i, j] for j in J) == 1, "Assign(%s)" % i)
for j in J:
model.addConstr(x[i, j] <= y[j], "Strong(%s,%s)" % (i, j))
model.addConstr(quicksum(y[j] for j in J) == k, "Facilities")
model.setObjective(quicksum(c[i, j] * x[i, j] for i in I for j in J), GRB.MINIMIZE)
model.update()
model.__data = x, y
return model
def distance(x1, y1, x2, y2):
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
def make_data(n, m, same=True):
"""
施設と顧客が同じ場合にはsameにTrueを，異なる場合にはFalseを入れる．
"""
if same == True:
I = range(n)
J = range(m)
x = [random.random() for i in range(max(m, n))]
y = [random.random() for i in range(max(m, n))]
else:
I = range(n)
J = range(n, n + m)
x = [random.random() for i in range(n + m)]
y = [random.random() for i in range(n + m)]
c = {}
for i in I:
for j in J:
c[i, j] = distance(x[i], y[i], x[j], y[j])
return I, J, c, x, y
random.seed(67)
n = 300
m = n
I, J, c, x_pos, y_pos = make_data(n, m, same=True)
k = 30
model = kmedian(I, J, c, k)
# model.Params.Threads = 1
model.Params.LogFile = "gurobi.log"
model.optimize()
EPS = 1.0e-6
x, y = model.__data
edges = [(i, j) for (i, j) in x if x[i, j].X > EPS]
facilities = [j for j in y if y[j].X > EPS]
print("Optimal value=", model.ObjVal)
print("Selected facilities:", facilities)
position = {}
for i in range(len(x_pos)):
position[i] = (x_pos[i], y_pos[i])
G = nx.Graph()
facilities = set(facilities)
unused = set(j for j in J if j not in facilities)
client = set(i for i in I if i not in facilities and i not in unused)
G.add_nodes_from(facilities)
G.add_nodes_from(unused)
for (i, j) in edges:
if i!=j:
G.add_edge(i, j)
nx.draw(G, position, with_labels=False, node_color="r", nodelist=facilities, node_size=80)
nx.draw(G, position, with_labels=False, node_color="c", nodelist=unused, node_size=30)
AMPLによる求解
AMPLモデル(kmedian.mod)
param n; # number or customers
param m; # number or facilities
param c {0..n-1, 0..m-1};
param k;
var x {0..n-1, 0..m-1} binary;
var y {0..m-1} binary;
minimize cost: sum {i in 0..n-1, j in 0..m-1} c[i,j] * x[i,j];
subject to
Service {i in 0..n-1}: sum {j in 0..m-1} x[i,j] = 1;
Kfacil: sum {j in 0..m-1} y[j] = k;
Activate {i in 0..n-1, j in 0..m-1}: x[i,j] <= y[j];ampl = AMPL(Environment(AMPL_PATH))
ampl.setOption('solver', 'highs')
ampl.read("../ampl/kmedian.mod")
ampl.param['n'] = n
ampl.param['m'] = n
ampl.param['k'] = k
ampl.param['c'] = c
ampl.solve()
cost = ampl.obj['cost']
print("Objective is:", cost.value())
facilities,edges = [],[]
x_ = ampl.var['x']
y_ = ampl.var['y']
容量制約付き施設配置問題
容量制約付き施設配置問題(capacitated facility location problem)は，以下のように定義される問題である．
顧客数を $n$， 施設数を $m$ とし， 顧客を $i=1,2,\ldots,n$， 施設を $j=1,2,\ldots,m$ と番号で表すものとする． また，顧客の集合を $I=\{1,2,\ldots,n\}$，施設の集合を $J=\{1,2,\ldots,m\}$ と記す．
顧客 $i$ の需要量を $d_i$，顧客 $i$ と施設 $j$ 間に $1$ 単位の需要が移動するときにかかる 輸送費用を $c_{ij}$，施設 $j$ を開設するときにかかる固定費用を $f_j$，容量を $M_j$ とする．
以下に定義される連続変数 $x_{ij}$ および $0$-$1$ 整数変数 $y_j$ を用いる．
$$ x_{ij}= 顧客 i の需要が施設 j によって満たされる量 $$$$ y_j = \left\{ \begin{array}{ll} 1 & 施設 j を開設するとき \\ 0 & それ以外のとき \end{array} \right. $$上の記号および変数を用いると，容量制約付き施設配置問題は以下の混合整数最適化問題 として定式化できる．
$$ \begin{array}{l l l } minimize & \sum_{j \in J} f_j y_j +\sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \sum_{i \in I} x_{ij} \leq M_j y_j & \forall j \in J \\ & x_{ij} \leq d_i y_j & \forall i \in I; j \in J \\ & x_{ij} \geq 0 & \forall i \in I; j \in J \\ & y_j \in \{ 0,1 \} & \forall j \in J \end{array} $$def flp(I, J, d, M, f, c):
"""flp -- model for the capacitated facility location problem
Parameters:
- I: set of customers
- J: set of facilities
- d[i]: demand for customer i
- M[j]: capacity of facility j
- f[j]: fixed cost for using a facility in point j
- c[i,j]: unit cost of servicing demand point i from facility j
Returns a model, ready to be solved.
"""
model = Model("flp")
x, y = {}, {}
for j in J:
y[j] = model.addVar(vtype="B", name="y(%s)" % j)
for i in I:
x[i, j] = model.addVar(vtype="C", name="x(%s,%s)" % (i, j))
model.update()
for i in I:
model.addConstr(quicksum(x[i, j] for j in J) == d[i], "Demand(%s)" % i)
for j in M:
model.addConstr(quicksum(x[i, j] for i in I) <= M[j] * y[j], "Capacity(%s)" % j)
for (i, j) in x:
model.addConstr(x[i, j] <= d[i] * y[j], "Strong(%s,%s)" % (i, j))
model.setObjective(
quicksum(f[j] * y[j] for j in J)
+ quicksum(c[i, j] * x[i, j] for i in I for j in J),
GRB.MINIMIZE,
)
model.update()
model.__data = x, y
return model
def make_data():
I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180}) # demand
J, M, f = multidict(
{1: [500, 1000], 2: [500, 1000], 3: [500, 1000]}
) # capacity, fixed costs
c = {
(1, 1): 4,
(1, 2): 6,
(1, 3): 9, # transportation costs
(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,
}
return I, J, d, M, f, c
I, J, d, M, f, c = make_data()
model = flp(I, J, d, M, f, c)
model.optimize()
EPS = 1.0e-6
x, y = model.__data
edges = [(i, j) for (i, j) in x if x[i, j].X > EPS]
facilities = [j for j in y if y[j].X > EPS]
print("Optimal value=", model.ObjVal)
print("Facilities at nodes:", facilities)
print("Edges:", edges)
非線形施設配置問題
容量制約付き施設配置問題において 倉庫における費用が出荷量の合計の凹費用関数になっている場合を考える．
ここでは施設 $j$ に対する固定費用 $f_j$ のかわりに， 施設 $j$ から輸送される量の合計 $X_j$ に対して以下の凹費用関数がかかるものと仮定する．
$$ f_j \sqrt{X_j} $$他の記号は，容量制約付き施設配置問題と同じである．
上の記号および変数を用いると，凹費用関数をもつ施設配置問題は以下の混合整数最適化問題 として定式化できる．
$$ \begin{array}{l l l } minimize & \sum_{j \in J} f_j \sqrt{\sum_{j \in J} x_{ij}} +\sum_{i \in I} \sum_{j \in J} c_{ij} x_{ij} & \\ s.t. & \sum_{j \in J} x_{ij} =d_i & \forall i \in I \\ & \sum_{i \in I} x_{ij} \leq M_j y_j & \forall j \in J \\ & x_{ij} \geq 0 & \forall i \in I; j \in J \end{array} $$まず，凸関数を最小化する問題を線形最適化問題に（近似的に）帰着する方法を考えよう． 凸関数を線形関数の繋げたものとして近似する．このように一部をみれば線形関数であるが， 繋ぎ目では折れ曲がった関数を区分的線形関数(piecewise linear function)と呼ぶ．
最初に，1変数の凸関数$f(x)$を区分的線形関数で近似する方法を考える．
変数 $x$ の範囲を $L \leq x \leq U$ とし， 範囲 $[L,U]$ を $K$ 個の区分 $[a_k,a_{k+1}] (k=0,1,\ldots,K-1)$ に分割する． ここで，$a_k$ は， $$ L= a_0 < a_1 < \cdots < a_{K-1} < a_K =U $$ を満たす点列である． 各点 $a_k$ における関数 $f$ の値を $b_k$ とする． $$ b_k =f(a_k) \ \ k=0,1,\ldots,K $$ 各区分 $[a_k,a_{k+1}]$ に対して，点 $(a_k,b_k)$ と点 $(a_{k+1},b_{k+1})$ を通る線分を引くことによって区分的線形関数を構成する．
線分は端点の凸結合で表現できる． $k$番目の点 $a_k$ に対する変数を$z_k$ とする． いま，$f(x)$ は凸関数であるので，隣り合う2つの添え字 $k, k+1$ に対してだけ $z_k$ が正となるので， $$ f(x) \approx \sum_{k=0}^K b_k z_k $$ $$ x= \sum_{k=0}^K a_k z_k $$ $$ \sum_{k=0}^K z_k=1 $$ $$ z_k \geq 0 \ \ \ \forall k=0,1,\ldots,K $$ が区分的線形近似を与える．
ここで考える平方根関数は凹関数なので，タイプ2の特殊順序集合(Special Ordered Set: SOS)と呼ばれる制約を付加するものである． タイプ2の特殊順序集合は，順序が付けられた集合（順序集合）に含まれる変数のうち， （与えられた順序のもとで）連続するたかだか2つが $0$ でない値をとることを規定する．
def convex_comb_sos(model, a, b):
"""convex_comb_sos -- add piecewise relation with gurobi's SOS constraints
Parameters:
- model: a model where to include the piecewise linear relation
- a[k]: x-coordinate of the k-th point in the piecewise linear relation
- b[k]: y-coordinate of the k-th point in the piecewise linear relation
Returns the model with the piecewise linear relation on added variables x, f, and z.
"""
K = len(a) - 1
z = {}
for k in range(K + 1):
z[k] = model.addVar(
lb=0, ub=1, vtype="C"
) # do not name variables for avoiding clash
x = model.addVar(lb=a[0], ub=a[K], vtype="C")
f = model.addVar(lb=-GRB.INFINITY, vtype="C")
model.update()
model.addConstr(x == quicksum(a[k] * z[k] for k in range(K + 1)))
model.addConstr(f == quicksum(b[k] * z[k] for k in range(K + 1)))
model.addConstr(quicksum(z[k] for k in range(K + 1)) == 1)
model.addSOS(GRB.SOS_TYPE2, [z[k] for k in range(K + 1)])
return x, f, z
def flp_nonlinear_sos(I, J, d, M, f, c, K):
"""flp_nonlinear_sos -- use model with SOS constraints
Parameters:
- I: set of customers
- J: set of facilities
- d[i]: demand for customer i
- M[j]: capacity of facility j
- f[j]: fixed cost for using a facility in point j
- c[i,j]: unit cost of servicing demand point i from facility j
- K: number of linear pieces for approximation of non-linear cost function
Returns a model, ready to be solved.
"""
a, b = {}, {}
for j in J:
U = M[j]
L = 0
width = U / float(K)
a[j] = [k * width for k in range(K + 1)]
b[j] = [f[j] * math.sqrt(value) for value in a[j]]
model = Model("nonlinear flp -- use model with SOS constraints")
x = {}
for j in J:
for i in I:
x[i, j] = model.addVar(
vtype="C", name=f"x(i,j)"
) # i's demand satisfied from j
model.update()
# total volume transported from plant j, corresponding (linearized) cost, selection variable:
X, F, z = {}, {}, {}
for j in J:
# add constraints for linking piecewise linear part:
X[j], F[j], z[j] = convex_comb_sos(model, a[j], b[j])
X[j].ub = M[j]
# constraints for customer's demand satisfaction
for i in I:
model.addConstr(quicksum(x[i, j] for j in J) == d[i], f"Demand(i)")
for j in J:
model.addConstr(quicksum(x[i, j] for i in I) == X[j], f"Capacity(j)")
model.setObjective(
quicksum(F[j] for j in J) + quicksum(c[i, j] * x[i, j] for j in J for i in I),
GRB.MINIMIZE,
)
model.update()
model.__data = x, X, F
return model
def distance(x1, y1, x2, y2):
return math.sqrt((x2 - x1) ** 2 + (y2 - y1) ** 2)
def make_data(n, m, same=True):
x, y = {}, {}
if same == True:
I = range(1, n + 1)
J = range(1, m + 1)
for i in range(1, 1 + max(m, n)): # positions of the points in the plane
x[i] = random.random()
y[i] = random.random()
else:
I = range(1, n + 1)
J = range(n + 1, n + m + 1)
for i in I: # positions of the points in the plane
x[i] = random.random()
y[i] = random.random()
for j in J: # positions of the points in the plane
x[j] = random.random()
y[j] = random.random()
f, c, d, M = {}, {}, {}, {}
total_demand = 0.0
for i in I:
for j in J:
c[i, j] = int(100 * distance(x[i], y[i], x[j], y[j])) + 1
d[i] = random.randint(1, 10)
total_demand += d[i]
total_cap = 0.0
r = {}
for j in J:
r[j] = random.randint(0, m)
f[j] = random.randint(100, 100 + r[j] * m)
M[j] = 1 + 100 + r[j] * m - f[j]
# M[j] = int(total_demand/m) + random.randint(1,m)
total_cap += M[j]
for j in J:
M[j] = int(M[j] * total_demand / total_cap + 1) + random.randint(0, r[j])
# print "%s\t%s\t%s" % (j,f[j],M[j])
# print "demand:",total_demand
# print "capacity:",sum([M[j] for j in J])
return I, J, d, M, f, c, x, y
def example():
I, d = multidict({1: 80, 2: 270, 3: 250, 4: 160, 5: 180}) # demand
J, M, f = multidict(
{10: [500, 100], 11: [500, 100], 12: [500, 100]}
) # capacity, fixed costs
c = {
(1, 10): 4,
(1, 11): 6,
(1, 12): 9, # transportation costs
(2, 10): 5,
(2, 11): 4,
(2, 12): 7,
(3, 10): 6,
(3, 11): 3,
(3, 12): 4,
(4, 10): 8,
(4, 11): 5,
(4, 12): 3,
(5, 10): 10,
(5, 11): 8,
(5, 12): 4,
}
x_pos = {
1: 0,
2: 0,
3: 0,
4: 0,
5: 0,
10: 2,
11: 2,
12: 2,
} # positions of the points in the plane
y_pos = {1: 2, 2: 1, 3: 0, 4: -1, 5: -2, 10: 1, 11: 0, 12: -1}
return I, J, d, M, f, c, x_pos, y_pos
I, J, d, M, f, c, x_pos, y_pos = example()
K = 4
model = flp_nonlinear_sos(I, J, d, M, f, c, K)
x, X, F = model.__data
model.Params.OutputFlag = 1 # silent/verbose mode
model.optimize()
edges = []
flow = {}
for (i, j) in sorted(x):
if x[i, j].X > EPS:
edges.append((i, j))
flow[(i, j)] = x[i, j].X
print("obj:", model.ObjVal, "\nedges", sorted(edges))
print("flows:", flow)
ハブ立地問題
ここでは，$p$-ハブ・メディアン問題 ($p$-hub median problem) に対する定式化を示す．
地点間には輸送量と距離（費用）が与えられており，各地点は，単一のハブに割り当てる必要があり，ハブ間は直送を行うものと仮定する． 目的は，総輸送費用の最小化である．
最初に記号を導入しておく．
- $n$: 点数
- $t_{ij}$: 地点 $i,j$ 間の輸送量
- $d_{ij}$: 地点 $i,j$ 間の距離
- $p$: 選択するハブの数
輸送費用は，距離に以下のパラメータを乗じて計算する．
- $\chi$: 地点からハブへの輸送
- $\alpha$: ハブ間の輸送
- $\delta$: ハブから地点への輸送
ベンチマーク問題例は，以下のサイトから入手できる．
http://people.brunel.ac.uk/~mastjjb/jeb/orlib/phubinfo.html
http://grafo.etsii.urjc.es/optsicom/uraphmp/
以下のデータ読み込みのコードでは，上のベンチマーク問題例が "../data/p-hub/" ディレクトリに保管されている場合を想定している． データの保管場所は適宜変更されたい．
folder = "../data/p-hub/"
f = open(folder + "AP40.txt")
data = f.readlines()
f.close()
n = int(data[0].split()[0])
t, d = {}, {}
data_ = []
for row in data[2:]:
data_.extend(list(map(float, row.split())))
count = 0
for i in range(n):
for j in range(n):
t[i, j] = data_[count]
count += 1
for i in range(n):
for j in range(n):
d[i, j] = data_[count]
count += 1
# 座標の読み込み
folder = "../data/p-hub/"
f = open(folder + "phub1.txt")
data2 = f.readlines()
f.close()
pos = {}
for i, row in enumerate(data2):
if row == "200\n":
break
for j in range(200):
xx, yy = list(map(int, data2[i + 1 + j].split()))
pos[j] = (xx, yy)
$p$-ハブ・メディアン問題に対する非凸2次関数定式化
最初の定式化は，非凸2次関数を目的関数としたものであるが，数理最適化ソルバーGurobiだと求解可能である．
以下の変数を用いる．
$x_{ik}$: 地点 $i$ がハブ $k$ に割り当てられるとき $1$ の $0$-$1$ 変数； ただし $x_{kk}=1$ のときは 地点 $k$ がハブであることを表す．
この変数を用いると，非凸2次な目的関数をもつ定式化は，以下のように書くことができる．
$$ \begin{array}{l l l} minimize & \sum_{i,j, k, \ell} t_{ij} ( \chi d_{ik} x_{ik} + \alpha d_{k \ell} x_{ik} x_{j\ell} +\delta d_{\ell j} x_{j\ell} ) & \\ s.t. & x_{ik} \leq x_{kk} & \forall i,k \\ & \sum_{k} x_{ik} = 1 & \forall i \\ & \sum_{k} x_{kk} = p & \\ & x_{ik} \in \{0,1\} & \forall i,k \end{array} $$p = 3
chi, alpha, delta = 3.0, 0.75, 2.0
N = list(range(n))
model = Model("p-hub")
x = {}
for i in N:
for k in N:
x[i, k] = model.addVar(vtype="B", name=f"x[{i},{k}]")
model.update()
for i in N:
model.addConstr(quicksum(x[i, k] for k in N) == 1)
model.addConstr(quicksum(x[k, k] for k in N) == p)
for i in N:
for k in N:
if k != i:
model.addConstr(x[i, k] <= x[k, k])
model.setObjective(
quicksum(sum(t[i, j] for j in N) * chi * d[i, k] * x[i, k] for i in N for k in N)
+ quicksum(
sum(t[i, j] for i in N) * delta * d[j, l] * x[j, l] for j in N for l in N
)
+ quicksum(
t[i, j] * alpha * d[k, l] * x[i, k] * x[j, l]
for i in N
for k in N
for j in N
for l in N
),
GRB.MINIMIZE
)
model.optimize()
G = nx.Graph()
for i, k in x:
if x[i, k].X > 0.001:
G.add_edge(i, k)
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=100, node_color="yellow")
plt.show()
$p$-ハブ・メディアン問題に対する線形定式化
ここでは，非凸2次関数に対応していない数理最適化ソルバーでも求解可能な定式化を示す．
以下の変数を用いる．
- $x_{ik}$:地点 $i$ がハブ $k$ に割り当てられるとき $1$ の $0$-$1$ 変数； ただし $z_{kk}=1$ のときは 地点 $k$ がハブであることを表す．
- $y_{ik\ell}$: 地点 $i$ で発生したものが，ハブ $k$ からハブ $\ell$ に輸送される量
最後の式は，地点 $i$ からハブ $k$ に運ばれた量に対するフロー保存式であり， これによって実数変数 $y_{ik\ell}$ が計算される．
p = 3
chi, alpha, delta = 3.0, 0.75, 2.0
N = list(range(n))
model = Model("p-hub")
x, y = {}, {}
for i in N:
for k in N:
x[i, k] = model.addVar(vtype="B", name=f"x[{i},{k}]")
for i in N:
for k in N:
for l in N:
y[i, k, l] = model.addVar(vtype="C", name=f"y[{i},{k},{l}]")
model.update()
for i in N:
model.addConstr(quicksum(x[i, k] for k in N) == 1)
model.addConstr(quicksum(x[k, k] for k in N) == p)
for i in N:
for k in N:
if k != i:
model.addConstr(x[i, k] <= x[k, k])
for i in N:
for k in N:
model.addConstr(
quicksum(y[i, k, l] for l in N) - quicksum(y[i, l, k] for l in N)
== (
sum(t[i, j] for j in N) * x[i, k]
- quicksum(t[i, j] * x[j, k] for j in N)
)
)
model.setObjective(
quicksum(sum(t[i, j] for j in N) * chi * d[i, k] * x[i, k] for i in N for k in N)
+ quicksum(sum(t[i, j] for i in N) * delta * d[l, j] * x[j, l] for j in N for l in N)
+ quicksum(alpha * d[k, l] * y[i, k, l] for i in N for k in N for l in N),
GRB.MINIMIZE
)
model.optimize()
G = nx.Graph()
for i, k in x:
if x[i, k].X > 0.001:
G.add_edge(i, k)
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=100, node_color="yellow")
plt.show()
$r$-割当 $p$-ハブ・メディアン問題
$r$-割当 $p$-ハブ・メディアン問題 ($r$-allocation $p$-hub median problem) とは， 各地点が高々 $r$ 個のハブに割当可能と仮定した拡張である． ．
まず，定式化に必要な記号を導入しておく．
- $n$: 点数
- $t_{ij}$: 地点 $i,j$ 間の輸送量
- $d_{ij}$: 地点 $i,j$ 間の距離
- $p$: 選択するハブの数
- $r$: 各点が接続可能なハブの数の上限
輸送費用は，距離に以下のパラメータを乗じて計算する．
- $\chi$: 地点からハブへの輸送
- $\alpha$: ハブ間の輸送
- $\delta$: ハブから地点への輸送
標準定式化
最初の定式化は変数は多いが，下界は強い標準的な定式化である．
以下の変数を用いる．
- $f_{ijk\ell}$: 地点 $i$ から地点 $j$ へ輸送する量が，パス $i \rightarrow k \rightarrow \ell \rightarrow j$ を通過する比率
- $z_{ik}$: 地点 $i$ がハブ $k$ に割り当てられるとき $1$ の $0$-$1$ 変数； ただし $z_{kk}=1$ のときは 地点 $k$ がハブであることを表す．
上の変数を使うと，$r$-割当 $p$-ハブ・メディアン問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i,j, k, \ell} t_{ij} ( \chi d_{ik}+ \alpha d_{k \ell} +\delta d_{\ell j} ) f_{ijk \ell} & \\ s.t. & \sum_{k} z_{ik} \leq r & \forall i \\ & z_{ik} \leq z_{kk} & \forall i,k \\ & \sum_{k} z_{kk} = p & \\ & \sum_{k, \ell} f_{ijk\ell} = 1 & \forall i,j \\ & \sum_{\ell} f_{ijk\ell} \leq z_{ik} & \forall i,j,k \\ & \sum_{k} f_{ijk\ell} \leq z_{j\ell} & \forall i,j,\ell \\ & f_{ijk \ell} \geq 0 & \forall i,j, k, \ell \\ & z_{ik} \in \{0,1\} & \forall i,k \end{array} $$p, r = 3, 2
chi, alpha, delta = 3.0, 0.75, 2.0
N = list(range(n))
model = Model("p-hub")
f, z = {}, {}
for i in N:
for j in N:
for k in N:
for l in N:
f[i, j, k, l] = model.addVar(vtype="C", name=f"f[{i},{j},{k},{l}]")
for i in N:
for k in N:
z[i, k] = model.addVar(vtype="B", name=f"z[{i},{k}]")
model.update()
for i in N:
model.addConstr(quicksum(z[i, k] for k in N) <= r)
for i in N:
for k in N:
if k != i:
model.addConstr(z[i, k] <= z[k, k])
model.addConstr(quicksum(z[k, k] for k in N) == p)
for i in N:
for j in N:
model.addConstr(quicksum(f[i, j, k, l] for k in N for l in N) == 1)
for i in N:
for j in N:
for k in N:
model.addConstr(quicksum(f[i, j, k, l] for l in N) <= z[i, k])
for i in N:
for j in N:
for l in N:
model.addConstr(quicksum(f[i, j, k, l] for k in N) <= z[j, l])
model.setObjective(
quicksum(
t[i, j] * (chi * d[i, k] + alpha * d[k, l] + delta * d[l, j]) * f[i, j, k, l]
for i, j, k, l in f
),
GRB.MINIMIZE,
)
model.optimize()
G = nx.Graph()
for i, k in z:
if z[i, k].X > 0.001:
G.add_edge(i, k)
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=100, node_color="yellow")
plt.show()
コンパクトな定式化
上の定式化は変数の数が $O(n^4)$ であり，強い下界を算出するが，大規模問題例には不向きである． 以下では $O(n^3)$ の変数を用いた定式化を示す．
以下の変数を用いる．
- $x_{ik}$: 地点 $i$ からハブ $k$ へ輸送する量
- $X_{i\ell j}$: 地点 $i$ から地点 $j$ への輸送が，ハブ $\ell$ から運ばれる量
- $y_{ik\ell}$: 地点 $i$ で発生したものが，ハブ $k$ からハブ $\ell$ に輸送される量
- $z_{ik}$: 地点 $i$ がハブ $k$ に割り当てられるとき $1$ の $0$-$1$ 変数； ただし $z_{kk}=1$ のときは 地点 $k$ がハブであることを表す．
上の変数を使うと，$r$-割当 $p$-ハブ・メディアン問題は以下のように定式化できる．
$$ \begin{array}{l l l} minimize & \sum_{i,k} \chi d_{ik} x_{ik} + \sum_{i,j,\ell} \delta d_{\ell j} X_{i\ell j}+ \sum_{i,k,\ell } \alpha d_{k \ell} y_{ik\ell } & \\ s.t. & \sum_{k} z_{ik} \leq r & \forall i \\ & z_{ik} \leq z_{kk} & \forall i,k \\ & \sum_{k} z_{kk} = p & \\ & \sum_{k} x_{ik} = \sum_{j} t_{ij} & \forall i \\ & x_{ik} \leq (\sum_{j} t_{ij}) z_{ik} & \forall i,k \\ & \sum_{\ell} X_{i\ell j} = t_{ij} & \forall i,j \\ & X_{i\ell j} \leq t_{ij} z_{j\ell } & \forall i,j,\ell \\ & \sum_{\ell} y_{ik\ell } -\sum_{\ell} y_{i\ell k} = x_{ik} - \sum_{j} X_{ikj} & \forall i,k \\ & x_{ik} \geq 0 & \forall i,k \\ & X_{i\ell j} \geq 0 & \forall i,\ell ,j \\ & y_{ik\ell } \geq 0 & \forall i,k,\ell \\ & z_{ik} \in \{0,1\} & \forall i,k \end{array} $$p, r = 3, 2
chi, alpha, delta = 3.0, 0.75, 2.0
N = list(range(n))
model = Model("p-hub")
x, X, y, z = {}, {}, {}, {}
for i in N:
for k in N:
x[i, k] = model.addVar(vtype="C", name=f"x[{i},{k}]")
for i in N:
for l in N:
for j in N:
X[i, l, j] = model.addVar(vtype="C", name=f"X[{i},{l},{j}]")
for i in N:
for k in N:
for l in N:
y[i, k, l] = model.addVar(vtype="C", name=f"y[{i},{k},{l}]")
for i in N:
for k in N:
z[i, k] = model.addVar(vtype="B", name=f"z[{i},{k}]")
model.update()
for i in N:
model.addConstr(quicksum(z[i, k] for k in N) <= r)
for i in N:
for k in N:
if k != i:
model.addConstr(z[i, k] <= z[k, k])
model.addConstr(quicksum(z[k, k] for k in N) == p)
for i in N:
model.addConstr(quicksum(x[i, k] for k in N) == sum(t[i, j] for j in N))
for i in N:
for k in N:
model.addConstr(x[i, k] <= sum(t[i, j] for j in N) * z[i, k])
for i in N:
for j in N:
model.addConstr(quicksum(X[i, l, j] for l in N) == t[i, j])
for i in N:
for j in N:
for l in N:
model.addConstr(X[i, l, j] <= t[i, j] * z[j, l])
for i in N:
for k in N:
model.addConstr(
quicksum(y[i, k, l] for l in N)
+ quicksum(X[i, k, j] for j in N)
- quicksum(y[i, l, k] for l in N)
- x[i, k]
== 0
)
model.setObjective(
quicksum(chi * d[i, k] * x[i, k] for i in N for k in N)
+ quicksum(delta * d[j, l] * X[i, l, j] for i in N for j in N for l in N)
+ quicksum(alpha * d[k, l] * y[i, k, l] for i in N for k in N for l in N),
GRB.MINIMIZE,
)
model.optimize()
G = nx.Graph()
for i, k in z:
if z[i, k].X > 0.001:
G.add_edge(i, k)
plt.figure()
nx.draw(G, pos=pos, with_labels=True, node_size=100, node_color="yellow")
plt.show()
ロジスティクス・ネットワーク設計問題
実際には，ロジスティクス・ネットワーク全体を最適化する必要がある．以下に簡単な例を示す．
集合・パラメータ・変数を以下に示す.
集合：
- $Cust$ : 顧客（群）の集合
- $Dc$ : 倉庫の集合
- $Plnt$ : 工場の集合
- $Prod$ : 製品（群）の集合
パラメータ：
需要量，固定費用，取扱量（生産量）上下限の数値は，単位期間（既定値は365日）に変換してあるものとする．
- $c_{ij}$ : 地点 $ij$ 間の1単位重量あたりの移動費用
- $w_{p}$ : 製品 $p$ の重量
- $d_{kp}$ : 顧客 $k$ における製品 $p$ の需要量
- $LB_j, UB_j$ : 倉庫 $j$ の取扱量の下限と上限；入庫する製品量の下限と上限を表す．
- $NLB, NUB$ : 開設する倉庫数の下限と上限
- $f_j$ : 倉庫 $j$ を開設する際の固定費用
- $v_j$ : 倉庫 $j$ における製品1単位あたりの変動費用
- $PLB_{ip}, PUB_{ip}$ : 工場 $i$ における製品 $p$ の生産量上限
変数：
$x_{ijp}$ : 地点 $ij$ 間の製品 $p$ のフロー量
$y_{j}$: 倉庫 $j$ を開設するとき $1$， それ以外のとき $0$ を表す$0$-$1$変数
定式化：
$$ \begin{array}{lll} minimize & \displaystyle\sum_{ i \in Plnt, j \in Dc, p \in Prod} (w_p c_{ij} + v_j) x_{ijp} + & \displaystyle\sum_{j \in Dc, k \in Cust, p \in Prod} w_p c_{jk} x_{jkp} + \displaystyle\sum_{j \in Dc} f_j y_j & \\ s.t. & \displaystyle\sum_{j \in Dc} x_{jkp} = d_{kp} & \forall k \in Cust, p \in Prod \\ & \displaystyle\sum_{ i \in Plnt} x_{ijp} = \sum_{k \in Cust} x_{jkp} & \forall j \in Dc, p \in Prod \\ & x_{jkp} \leq d_{kp} y_j & \forall j \in Dc, k \in Cust, p \in Prod \\ & LB_j y_j \leq \displaystyle\sum_{i \in Plnt, p \in Prod} x_{ijp} \leq UB_j y_j & \forall j \in Dc \\ & \displaystyle\sum_{j \in Dc} x_{ijp} \leq PUB_{ip} & \forall i \in Plnt, p \in Prod \\ & NUB \leq \displaystyle\sum_{j \in Dc} y_j \leq NUB & \end{array} $$上のモデルの拡張を組み込んだロジスティクス・ネットワーク最適化システムとしてMELOS (MEta Logistic network Optimization System: https://www.logopt.com/melos/ ) が開発されている．