おまじないとして以下のコードをコピペして実行してください.
def draw_region(_demands, _pfls):
for i in range(len(_demands)):
plt.text(x=_demands[i,0]-0.13, y=_demands[i,1]-0.15, s=str(i+1), fontdict=dict(color='white',size=10))
plt.text(x=_demands[i,0]+0.1, y=_demands[i,1]+0.1, s=str(int(_demands[i,2])), fontdict=dict(color='blue',size=10))
plt.scatter(_demands[:,0], _demands[:,1], s=150)
plt.scatter(_pfls[:,0], _pfls[:,1], c='red', s=200)
for i in range(len(_pfls)):
plt.text(x=_pfls[i,0]-0.12, y=_pfls[i,1]-0.15, s=str(i+1), fontdict=dict(color='black',size=10))
plt.show()
def draw_facility(_demands, _pfls, _y):
for i in range(len(_demands)):
plt.text(x=_demands[i,0]-0.13, y=_demands[i,1]-0.15, s=str(i+1), fontdict=dict(color='white',size=10))
plt.text(x=_demands[i,0]+0.1, y=_demands[i,1]+0.1, s=str(int(_demands[i,2])), fontdict=dict(color='blue',size=10))
plt.scatter(_demands[:,0], _demands[:,1], s=150)
plt.scatter(_pfls[:,0], _pfls[:,1], c='red', s=200)
x_coord = []
y_coord = []
for i, val in enumerate(_y):
if val == 1:
x_coord.append(_pfls[i,0])
y_coord.append(_pfls[i,1])
plt.scatter(x_coord, y_coord, c='yellow', s=200)
for i in range(len(_pfls)):
plt.text(x=_pfls[i,0]-0.12, y=_pfls[i,1]-0.15, s=str(i+1), fontdict=dict(color='black',size=10))
def results(_demands, _pfls, _x, _y):
x_ori = []
x_des = []
y_ori = []
y_des = []
for i, row in enumerate(_x):
for j, val in enumerate(row):
if val == 1:
x_ori.append(demands[i,:2][0])
x_des.append(pfls[j][0])
y_ori.append(demands[i,:2][1])
y_des.append(pfls[j][1])
plt.plot([x_ori, x_des], [y_ori, y_des], color='blue')
for i in range(len(demands)):
plt.text(x=_demands[i,0]-0.13, y=_demands[i,1]-0.15, s=str(i+1), fontdict=dict(color='white',size=10))
plt.text(x=_demands[i,0]+0.1, y=_demands[i,1]+0.1, s=str(int(_demands[i,2])), fontdict=dict(color='blue',size=10))
plt.scatter(_demands[:,0], _demands[:,1], s=150)
x_coord = []
y_coord = []
for i, val in enumerate(_y):
if val == 1:
x_coord.append(_pfls[i,0])
y_coord.append(_pfls[i,1])
plt.scatter(_pfls[:,0], _pfls[:,1], c='red', s=200)
plt.scatter(x_coord, y_coord, c='yellow', s=200)
for i in range(len(_pfls)):
plt.text(x=_pfls[i,0]-0.12, y=_pfls[i,1]-0.15, s=str(i+1), fontdict=dict(color='black',size=10))
# google collabを利用している人は、はじめに以下を実行してください。
!pip install pulp
12個の街と3つの店舗配置候補点が存在する地域を考えてみましょう. まずは街と店舗配置候補点の座標が書き込まれているデータをロードしましょう.
import numpy as np
!wget http://int-info.com/PyLearn/Grad_4/data_for_study/pfl.csv
!wget http://int-info.com/PyLearn/Grad_4/data_for_study/demand.csv
# 街の座標
demands = np.loadtxt('demand.csv', delimiter=',')
# 店舗配置候補点の座標
pfls = np.loadtxt('pfl.csv', delimiter=',')
# 街から店舗配置候補点までの距離
distance = np.zeros([len(demands), len(pfls)])
import matplotlib.pyplot as plt
draw_region(demands, pfls)
この画像におけるそれぞれの点をノードと呼びます. 青いノードが街,赤いノードが店舗配置候補点です. 今後,街のことは需要点(店の需要である住民がいるため),店舗配置候補点を配置点と呼びます. 需要点と配置点の中に書かれている数字は,そのノードの番号を示しています. また,需要点の右上に書かれている数字はその需要点の人口を示しています.
ここでこの地域にスーパーマーケットをいくつか出店することを考えてみましょう. 例えば,2つのスーパーマーケットを出店する場合にはどこに店舗を配置すれば効率がよいのでしょうか. この効率を考えるためにあらかじめ各需要点から各配置点までの距離を計算しておきましょう. 以下のコードで計算できます.
def cal_distance(_demands, _pfls):
_distance = np.zeros([len(_demands), len(_pfls)])
for i, demand in enumerate(_demands):
for j, pfl in enumerate(_pfls):
p1 = pfl
p2 = demand[:2].astype(int)
_distance[i][j] = np.sqrt((p1[0]-p2[0])**2 + (p1[1]-p2[1])**2)
return _distance
distance = cal_distance(demands, pfls)
p-median問題 とはざっくりいうと,対象地域の全住民の移動距離が最も小さくなるような p個の施設を配置する問題,です. またp-median問題は目的関数と制約条件から構成されます.
目的関数とは最小化もしくは最大化すべき指標のことで,制約条件とは目的関数を最小化もしくは最大化する際に満たすべき条件のことです. ここを抑えて p-median 問題の目的関数と制約条件を見てみましょう.
目的関数(最小化)
制約条件
このように目的関数と制約条件で表された問題を最適化問題と呼びます.
重み付き総移動距離とは街の住民とその住民が利用する店舗までの距離の総和です. これは"街の人口×その街から店舗が配置された街までの距離" で表すことができます. 例えば街Aの人口が100人,街Bの人口が200人であり,店舗までの距離がそれぞれ 1km と2km であれば
Aの人口×Aから店舗までの距離 + Bの人口×Bから店舗までの距離
が重み付き移動距離となります.
前項では文章でいろいろと説明してきましたが,コンピュータで扱えるように数学でこの問題を記述してみましょう.
まず,二つの記号を用意します. 一つは店舗を配置することのできるノードを表す$i$ です(これを配置点と呼びます). もう一つは住民が住んでいるノードを表す$j$ です(これを需要点と呼びます). 今回のグラフでは$i$ が3個,$j$ が12個ですから
$i = 1, 2, \cdots, 12$
$j = 1, 2, 3$
となります.
次に需要点$j$ の住民が配置点$i$ に配分されているかどうかを表す変数$x_{ij}$ を考えます. この変数$x_{ij}$ は0 か1 のみを取る変数であり,例えば,需要点3 の住民が配置点2の店舗を利用するとき$x_{3, 2}=1$ となります. もう少し一般的に書いてみると次のようになります.
$x_{ij} = 1:$ 需要点$i$ の住民が配置点$j$ の店舗を利用する,$0$: その他
次に配置点$i$ に施設が配置されているかどうかを表す変数$y_i$ を考えます. 変数$y_i$ も同様に0 か1 のみを取る変数となっています. すなわち,
$y_j = 1:$ 配置点$j$ に店舗が配置されている,$0$: その他
となります.
あらかじめ与えられた値は定数として扱います. 今回の例ではそれぞれの街の人口と各街同士の距離が定数になります. それぞれ次のように定義しましょう.
目的関数は前項で述べたように重み付き総移動距離,すなわち人口×移動距離の総和です. これを上述した記号や変数,定数で表現してみましょう. 需要点$i$ から配置点$j$ への重み付き総移動距離は次のように表現できます.
$p_i \cdot d_{ij}$
この式は,需要点$i$ の住民が配置点$j$ に配置された店舗を利用する場合,すなわち$x_{ij} = 1$ のときのみ有効にならなければなりません.すなわち
$p_i \cdot d_{ij} \cdot x_{ij}$
となります. この式が需要点$i$ から配置点$j$ に対する重み付き総移動距離となります. 今度は需要点$i$ からすべての配置点に対する重み付き総移動距離を考えてみましょう. 上記の式をすべての配置点について足してあげれば実現できるはずです. すなわち,
$p_i \cdot \sum_{j = 1}^{3} d_{ij} \cdot x_{ij}$
となります. 簡単な復習になりますが$\sum$ は総和を意味しており,例えば $\sum_{i = 1}^{3} a_i = a_1 + a_2 + a_3$ です. これで街$i$ から配置点に対する重み付き総移動距離を数式で表現sるうことができました. 最後に需要点$i$ だけでなくすべての需要点について同じような重み付き総移動距離を考えてそれらすべてを足してみましょう. すなわち,
$\sum_{i = 1}^{12} p_i \cdot \sum_{j = 1}^{3} d_{ij} \cdot x_{ij}$
です. これが対象地域上の全住民の移動距離を考慮した重み付き総移動距離となり,目的関数となるわけです. 今回は目的関数を最小化するので,これを以下のように表現します.
$\sum_{i = 1}^{12} p_i \cdot \sum_{j = 1}^{3} d_{ij} \cdot x_{ij} \rightarrow \min$
一見難しい表記に見えますが$\rightarrow \min$ は左の目的関数を最小化するという意味です.
今回の問題では制約条件が三つありました.それぞれ順番に数学的に表現してみましょう.
まずは「各需要点の住民が利用する店舗は一つ」です. 住民が店舗を利用するかどうかは変数$x_{ij}$ で表現することができました. そこで変数$x_{ij}$ を使ってこの制約条件を表現してみましょう. まずは具体例として「需要点$1$ の住民が利用する店舗は一つだけ」という制約条件を考えましょう. これは
$\sum_{j=1}^{3} x_{1,j} = 1$
となります. $x_{1,j} = 1$ は需要点$1$ の住民が配置点$j$ の店舗を利用する,という意味なのですべての配置店$j$ について$x_{1,j}$ を足してあげればいいわけです. この式は需要点$1$ についてのみ考えているので,すべての需要点について考えてあげると以下のように9個の式ができます.
$\sum_{j=1}^{3} x_{1,j} = 1$
$\sum_{j=1}^{3} x_{2,j} = 1$
$\cdots$
$\sum_{j=1}^{3} x_{12,j} = 1$
これらはまとめて次のように表記できます.
$\sum_{j=1}^{3} x_{i,j} = 1, \forall i\in \{1, 2, \cdots, 12\}$
$\forall i\in \{1, 2, \cdots, 12\}$ はこういう表記なんだと納得してください.
次に「配置点$j$ に店舗が配置されていなければ需要点の住民は配置点$j$に買い物に行けない」です. 当たり前のことですが,これもしっかりと数式で表してあげます. 繰り返しになりますが,住民が店舗を利用するかどうかは変数$x_{ij}$ で表せます. そのため,$y_{j}=1$ (店舗が配置される)でなければ$x_{ij}=1$ (その店舗を利用する)になりません. つまり制約条件は次のような式で表せます.
$x_{ij} \leq y_{j}$
配置点$j$ に店舗が配置された($y_{j}=1$)時のみ需要点$i$ の住民がその店舗を利用できる$x_{ij}=1$,という意味になっています. この式をすべての需要点と配置点について考えると次のように表記できます.
$x_{ij} \leq y_{j}, \forall i \in \{1, 2, \cdots, 12\}, \forall j \in \{1, 2, 3\}$
これはこのような表記なんだと納得してください.
次に「出店されるスーパーマーケットの数は2店舗」です. スーパーマーケットが配置店$j$ に出店されているかどうかは変数$y_{j}$ で表現することができました. 変数$y_{j}$ が$1$ であれば出店されており,$0$ であれば出店されていない,という意味になります. そのため,出店されているスーパーマーケットの数は変数$y_{j}$ の足し算で考えることができます. すなわち出店されているスーパーマーケットの数は
$\sum_{j=1}^{3} y_{j}$
と表現することができます. 今回は出店舗数が$2$ ですからこの制約条件は以下のように表現できます.
$\sum_{j=1}^{3} y_{j} = 2$
これで出店舗数に関する制約条件を数学的に表現することができました.
目的関数と制約条件を数学的に表現することができたので,これらをまとめて数理モデルとして表現しましょう. p-median問題の数理モデルは次のように表現されます.
\begin{align} &\sum_{i = 1}^{12} p_i \cdot \sum_{j = 1}^{3} d_{ij} \cdot x_{ij} \rightarrow \min, \label{eq:Obj} \\ &\text{s.t. } \ \ \sum_{j=1}^{3} x_{i,j} = 1, \forall i\in \{1, 2, \cdots, 12\}, \\ &\ \ \ \ \ \ \ \ x_{ij} \leq y_{j}, \forall i \in \{1, 2, \cdots, 12\}, \forall j \in \{1, 2, 3\}, \\ &\ \ \ \ \ \ \ \ \sum_{j=1}^{3} y_{j} = 2. \end{align}それではここからはこれをコンピュータに解かせてみましょう.
まずは需要点と配置点の数や距離,需要点の人口などの定数を定義していきます.
# 需要点の数
demand = len(demands)
# 配置点の数
pfl = len(pfls)
# 各需要点の人口
population = demands[:,2]
# 各需要点と各配置点の距離
distance = cal_distance(demands, pfls)
# 配置する施設の数
stores = 2
各需要点の人口を確認してみましょう.
list(population)
一番上の行が需要点1,二つ目が2, といった感じになっています. つまり需要点1 の人口は500人,需要点2 は300人となっていて,最初に示した図と一致していることがわかります.
次に各需要点と各配置点の距離を確認してみましょう.
distance
一番目の行の一列目は需要点1 から施設点1 への距離を表しています. 距離の単位は何でもよいですが,ここでは[km] としておきましょう. つまり需要点1 から配置点1 への距離は約7.6[km] であるということがわかります.
これらの定数をもとに今回は'pulp' というソフトを用いてp-median問題を定式化して解を求めてみます. 以下のコードが先ほど示した数理モデルを表しています. 余裕のある人はコードの意味を自分で調べてみるといいかもしれません.
import pulp
def solve_pmedian(_demand, _pfl, _population, _distance, _stores):
solver = pulp.PULP_CBC_CMD(msg=True, logPath='')
# 問題設定
slp = pulp.LpProblem('p-median', pulp.LpMinimize)
# 変数設定
x = [[pulp.LpVariable('x%d %d'%(i+1, j+1), cat=pulp.LpBinary) for j in range(_pfl)] for i in range(_demand)]
y = [pulp.LpVariable('y%d'%(j+1), cat=pulp.LpBinary) for j in range(_pfl)]
# 目的関数
slp += pulp.lpSum([_population[i] * pulp.lpSum([_distance[i, j] * x[i][j] for j in range(_pfl)]) for i in range(_demand)])
# 制約条件
for i in range(_demand):
# 制約条件1
slp += pulp.lpSum([x[i][j] for j in range(pfl)]) == 1
for j in range(_pfl):
# 制約条件2
slp += x[i][j] - y[j] <= 0
# 制約条件3
slp += pulp.lpSum([y[j] for j in range(_pfl)]) == _stores
# 求解
slp.solve(solver)
# 目的関数値保存
_f = pulp.value(slp.objective)
# 変数x とy を保存
_x = np.array([[pulp.value(x[i][j]) for j in range(_pfl)] for i in range(_demand)], dtype=np.int)
_y = np.array([pulp.value(y[j]) for j in range(_pfl)], dtype=np.int)
return _f, _x, _y
f, x, y = solve_pmedian(demand, pfl, population, distance, stores)
'Result - Optimal solution found' という行があります.
これは最適解,すなわち最も目的関数が小さくなる解を見つけることができた,ということを示しています.
それでは,実際にどこに店舗が配置されたのか確認してみましょう.
店舗の配置を示す変数はy
でしたのでy
の中身を確認してみます.
y
これは$y_1 = 1, y_2 = 1, y_3 = 0$ であることを示しています. すなわち,配置点1 と2 に店舗が配置されたということになります.
次に実際に重み付き総移動距離はいくつくらいになったのか確認してみましょう.
print(f'{f:.2f}')
これは対象地域上の住民が今回配置された店舗を使ったとき,その総移動距離が18937.25[km] になることを表しています.
それでは図に示して確認してみましょう. 施設が配置されたノードを黄色で示してみます.
import matplotlib.pyplot as plt
draw_facility(demands, pfls, y)
pfls[:2,0]
配置されたノードを見てみると周りに需要点がきれいに散らばっており,なんとなく正しそうですね.
もっと具体的にどの需要点がどの店舗を利用しているのか図示して確認してみましょう. これには変数$x_{ij}$ の情報を使う必要があります. まずは中身を見てみましょう.
x
例えば1行目の2列目が1
になっていますが,これは需要点1 の住民が配置点2 に配置された店舗を利用していることを示しています.
この情報を使えばだれがどの店舗を使っているのか図示できそうですね.
import numpy as np
import matplotlib.pyplot as plt
results(demands, pfls, x, y)
詳細に確認することは難しいですが,配置点から需要点までの距離のバランスがよくうまく配置されていることがわかります.
ちなみに重み付き総移動距離とは「需要点と配置点を結ぶ線の長さ×その需要点の人口」をすべての線について足したものです. 図で見ると少しわかりやすいでしょうか.
(1) 今回はあらかじめ需要点と配置点が用意されていました.課題では配置点の数と座標を変えてみて新たな地域を作成してください.
座標はx座標,y座標とも1~10です. 例えば,座標(1, 3),(5, 5)に新たな配置点を置きたければ次のようにして下さい.
pfls = [] # この行はいじらない
pfls.append([1, 3])
pfls.append([5, 5])
# ↑に pfls.append([座標])と入力することで点を追加できる.
pfls = np.array(pfls) # この行はいじらない
新たに設定した座標を確認したければ以下のようにコマンドを打ってください.
draw_region(demands, pfls)
(2) 新たに作成した地域で問題を作成するために必要な定数を定義してください.なお,配置する店舗の数は自由に設定してください.
(3) 配置する店舗の数が変わっているため数理モデルの表記も変わります.このページ内の「p-median 問題の数理モデル」を参考に作成した問題を表す数理モデルを表記してください.
(4) 新たに定義した定数を使って問題を解いてみてください.
問題は次のコードで解くことができ,目的関数値と変数$x_{ij}, y_j$ はそれぞれf
,x
,y
に格納されます.
f, x, y = solve_pmedian(demand, pfl, population, distance, stores)
(5) 得られた結果を図示してください.
次のコードで図示できます.
results(demands, pfls, x, y)
参考にした資料を、2〜3件記載してください。以下、書き方の例です。
本文中で引用する場合は「大前らはXXXを実施している [1]。また、〜」のように、どこで引用したのか明白にしてください。