import random
import math
from gurobipy import Model, quicksum, GRB
import networkx as nx
from scipy.optimize import minimize
from numpy import array
Weber問題
平面上のどこにでも施設を配置できるタイプの立地問題は，Weber問題(Weber problem)とよばれる．
各顧客 $i$ は平面上に分布しているものとし，その座標を $(x_i,y_i)$ とする． 顧客は需要量 $w_i$ をもち，目的関数は施設と顧客の間の距離に 需要量を乗じたものの和とする． 顧客と施設間の距離は直線距離とする．
顧客の集合 $I$ に対して，単一の施設の配置地点 $(X,Y)$ を決定する問題は， $$ f(X,Y) =\sum_{i \in I} w_i \sqrt{ (X-x_i)^2 +(Y-y_i)^2 } $$ を最小にする $(X,Y)$ を求める問題になる． これは，不動点アルゴリズムを用いることによって容易に求解できるが，ここでは数理最適化ソルバーで解いてみよう．
顧客 $i$ と施設の間の直線距離を表す変数 $z_i$ を導入すると， Weber問題は次のように定式化することができる．
$$ \begin{array}{ll} minimize & \sum_{i\in I} w_i z_i & \\ s.t. & (x_i-X)^2 + (y_i-Y)^2 \leq z_i^2 & \forall i \in I \\ \end{array} $$上の問題の制約は，2次錐制約(second-order cone constraint)と呼ばれ，モダンな数理最適化ソルバー（たとえばGurobi）を使うと効率良く処理することができる．
def weber(I, x, y, w):
"""weber: model for solving the single source weber problem using soco.
Parameters:
- I: set of customers
- x[i]: x position of customer i
- y[i]: y position of customer i
- w[i]: weight of customer i
Returns a model, ready to be solved.
"""
model = Model("weber")
X, Y, z, xaux, yaux = {}, {}, {}, {}, {}
X = model.addVar(lb=-GRB.INFINITY, vtype="C", name="X")
Y = model.addVar(lb=-GRB.INFINITY, vtype="C", name="Y")
for i in I:
z[i] = model.addVar(vtype="C", name="z(%s)" % (i))
xaux[i] = model.addVar(lb=-GRB.INFINITY, vtype="C", name="xaux(%s)" % (i))
yaux[i] = model.addVar(lb=-GRB.INFINITY, vtype="C", name="yaux(%s)" % (i))
model.update()
for i in I:
model.addConstr(
xaux[i] * xaux[i] + yaux[i] * yaux[i] <= z[i] * z[i], "MinDist(%s)" % (i)
)
model.addConstr(xaux[i] == (x[i] - X), "xAux(%s)" % (i))
model.addConstr(yaux[i] == (y[i] - Y), "yAux(%s)" % (i))
model.setObjective(quicksum(w[i] * z[i] for i in I), GRB.MINIMIZE)
model.update()
model.__data = X, Y, z
return model
def make_data(n, m):
I = range(1, n + 1)
J = range(1, m + 1)
x, y, w = {}, {}, {}
for i in I:
x[i] = random.randint(0, 100)
y[i] = random.randint(0, 100)
w[i] = random.randint(1, 5)
return I, J, x, y, w
random.seed(3)
n = 7
m = 1
I, J, x, y, w = make_data(n, m)
model = weber(I, x, y, w)
model.optimize()
X, Y, z = model.__data
G = nx.Graph()
G.add_nodes_from(I)
G.add_nodes_from(["D"])
position = {}
for i in I:
position[i] = (x[i], y[i])
position["D"] = (round(X.X), round(Y.X))
nx.draw(G, pos=position, node_size=200, node_color="g", nodelist=I)
nx.draw(G, pos=position, node_size=400, node_color="red", nodelist=["D"], alpha=0.5)
SciPyを用いて最適化
SciPyの最適化モジュール optimize の minimize(f,x0) を用いて非線形最適化を行うこともできる． ここでfは最小化する関数であり，x0は初期点である．
引数methodで探索のためのアルゴリズムを設定することができる．以下のアルゴリズムが準備されている． ここでは 'SLSQP'を用いる．
'Nelder-Mead'：Nelder–Mead法(単体法)
'Powell'：Powellの共役方向法
'CG'：共役勾配法(conjugate gradient)
'BFGS'：Broyden–Fletcher–Goldfarb–Shanno(BFGS)法
'Newton-CG'：Newton共役勾配法
'L-BFGS-B'：記憶制限付きBFGS法
'TNC'：打ち切りNewton共約勾配法
'COBYLA'：Constrained Optimization BY Linear Approximation
'SLSQP'：Sequential Least SQuares Programming
'dogleg'：ドッグレッグ信頼領域法
'trust-ncg'：信頼領域Newton共約勾配法
random.seed(3)
n = 7
m = 1
I, J, x, y, w = make_data(n, m)
X = array(1)
def f(X):
total_cost = 0.0
for i in I:
total_cost += w[i] * math.sqrt(((x[i] - X[0]) ** 2 + (y[i] - X[1]) ** 2))
return total_cost
res = minimize(f, (80, 10), method="SLSQP", constraints=())
print("Result=",res)
G = nx.Graph()
G.add_nodes_from(I)
G.add_nodes_from(
"%s" % j for j in J
)
position = {}
for i in I:
position[i] = (x[i], y[i])
for j in J:
position["%s" % j] = (res.x[0], res.x[1])
nx.draw(G, position, node_color="g", nodelist=I, with_labels=False)
nx.draw(G, position, node_color="r", nodelist=["%s" % j for j in J], alpha=0.5)
複数施設Weber問題
Weber問題の拡張として，複数の施設を平面上に立地する場合を考える．
各顧客は，複数ある施設のうち，最も近い地点に割り振られるものとする． 施設の集合 $J$ を与えたとき， 各施設 $j (\in J$の配置地点 $(X_j,Y_j))$ を決定する問題は， $$ f(X,Y) =\sum_{i \in I} w_i \min_{j \in J} \sqrt{ (X_j-x_i)^2 +(Y_j-y_i)^2 } $$ を最小にする $(X_j,Y_j)$ を求める問題になる．
この問題の変数は，施設の座標 $(X_j, Y_j) (j \in J)$ の他に，以下の変数を用いる．
- $v_{ij}:$ 顧客 $i$ と施設 $j$ の間の直線距離
- $u_{ij}:$ 顧客 $i$ が施設 $j$ に割り当てられるとき $1$ になる $0$-$1$ 変数
これらの変数を用いると， 複数施設Weber問題は次のように定式化することができる．
$$ \begin{array}{ll} minimize & \sum_{i\in I} w_i u_{ij} v_{ij} & \\ s.t. & (x_i-X_j)^2 + (y_i-Y_j)^2 \leq v_{ij}^2 & \forall i \in I, j \in J \\ & \sum_{j \in J} u_{ij} = 1 & \forall i \in I \\ & u_{ij} \in \{0,1 \} & \forall i \in I, j \in J \end{array} $$上の問題の目的関数は，非凸の2次関数であるが，モダンな数理最適化ソルバー（たとえばGurobi）だと，自動的に変形して求解してくれる．
ただし，大規模問題例に対しては不向きである．実際には，Weber問題を解くための反復解法 (Weiszfeld法） を用いて，近似的に最適な施設を求める方法がある． それらを組み込んだ連続型施設最適化のためのシステムとしてMELOS-GF (MEta Logistic network Optimization System Green Field: https://www.logopt.com/melos/ ) が開発されている．
def weber_MS(I, J, x, y, w):
"""weber -- model for solving the weber problem using soco (multiple source version).
Parameters:
- I: set of customers
- J: set of potential facilities
- x[i]: x position of customer i
- y[i]: y position of customer i
- w[i]: weight of customer i
Returns a model, ready to be solved.
"""
M = max([((x[i] - x[j]) ** 2 + (y[i] - y[j]) ** 2) for i in I for j in I])
model = Model("weber - multiple source")
X, Y, v, u = {}, {}, {}, {}
xaux, yaux, uaux = {}, {}, {}
for j in J:
X[j] = model.addVar(lb=-GRB.INFINITY, vtype="C", name="X(%s)" % j)
Y[j] = model.addVar(lb=-GRB.INFINITY, vtype="C", name="Y(%s)" % j)
for i in I:
v[i, j] = model.addVar(vtype="C", name="v(%s,%s)" % (i, j))
u[i, j] = model.addVar(vtype="B", name="u(%s,%s)" % (i, j))
xaux[i, j] = model.addVar(
lb=-GRB.INFINITY, vtype="C", name="xaux(%s,%s)" % (i, j)
)
yaux[i, j] = model.addVar(
lb=-GRB.INFINITY, vtype="C", name="yaux(%s,%s)" % (i, j)
)
model.update()
for i in I:
model.addConstr(quicksum(u[i, j] for j in J) == 1, "Assign(%s)" % i)
for j in J:
model.addConstr(
xaux[i, j] * xaux[i, j] + yaux[i, j] * yaux[i, j] <= v[i, j] * v[i, j],
"MinDist(%s,%s)" % (i, j),
)
model.addConstr(xaux[i, j] == (x[i] - X[j]), "xAux(%s,%s)" % (i, j))
model.addConstr(yaux[i, j] == (y[i] - Y[j]), "yAux(%s,%s)" % (i, j))
model.setObjective(
quicksum(w[i] * v[i, j] * u[i, j] for i in I for j in J), GRB.MINIMIZE
)
model.update()
model.__data = X, Y, v, u
return model
random.seed(3)
n = 15
m = 3
I, J, x, y, w = make_data(n, m)
model = weber_MS(I, J, x, y, w)
model.optimize()
X, Y, v, u = model.__data
G = nx.Graph()
G.add_nodes_from(I)
J_list = []
for j in J:
G.add_node(f"D{j}")
J_list.append(f"D{j}")
for i in I:
for j, n in enumerate(J_list):
if u[i, j + 1].X > 0.001:
G.add_edge(i, n)
position = {}
for i in I:
position[i] = (x[i], y[i])
for j, n in enumerate(J_list):
position[n] = (X[j + 1].X, Y[j + 1].X)
nx.draw(G, pos=position, node_size=200, with_labels=True, node_color="y", nodelist=I)
nx.draw(
G,
pos=position,
node_size=400,
with_labels=True,
node_color="red",
nodelist=J_list,
alpha=0.5,
)
