これなに
4月から新しく下宿して、通勤・通学を始めた方も多いかもしれません。
東京だと、なかなか近くに住むのも大変です。
さて、みんなの住むところをシャッフルして、自由に引っ越し可能とした場合、総通勤時間を最小化するとどうなるでしょう？
(勤務地も住宅の位置も同じで、住む人を変えるものとします)
解いてみる
通勤時間の分布ですが、通勤に関するアンケート調査結果の電車通勤を見ると、平均 33.23分、標準偏差 21.80分とありますので、それを使います。分布は、適当に対数正規分布にしましょう。
住宅は、半径60分の円内に一様にランダムに作ります。勤務地は、そこから方向はランダムで、上記の距離離れたところとします。
Pythonで表示してみましょう。
python
%matplotlib inline
import numpy as np, matplotlib.pyplot as plt, networkx as nx
from pulp import *
from ortoolpy import addbinvars
plt.rcParams['font.family'] = 'IPAexGothic'
n = 100 # 100人分
np.random.seed(1)
dirc1 = np.random.uniform(0,np.pi*2,n)
dist1 = np.random.triangular(0,60,60,n)
dist2 = np.random.lognormal(3.2644,0.5838,n)
dirc2 = np.random.uniform(0,np.pi*2,n)
xs = np.sin(dirc1)*dist1 # 住宅X座標
ys = np.cos(dirc1)*dist1 # 住宅y座標
xt = np.sin(dirc2)*dist2 + xs # 勤務地X座標
yt = np.cos(dirc2)*dist2 + ys # 勤務地y座標
print('%.2f %.2f'%(dist2.mean(), dist2.std()))
plt.title('元の通勤時間分布')
plt.hist(a,bins=100, range=(0,120));
>>>
33.23 21.80
python
plt.title('住宅')
plt.scatter(xs,ys);
python
plt.title('勤務地')
plt.scatter(xt,yt);
python
print(dist2.sum()) # 元の総移動時間
>>>
3323.3549788106984
最適化してみます。
python
g = nx.Graph()
for i1,(x1,y1) in enumerate(zip(xs,ys)):
for i2,(x2,y2) in enumerate(zip(xt,yt)):
g.add_edge(i1, i2+n, weight=1000-np.linalg.norm([x1-x2,y1-y2]))
m = LpProblem(sense=LpMaximize)
x = np.array(addbinvars(n,n))
m += lpDot(x.flatten(),[g.edge[i][j]['weight'] for i,j in g.edges()])
for i in range(n):
m += lpSum(x[i]) == 1
m += lpSum(x[:,i]) == 1
%time m.solve()
s = (np.vectorize(value)(x)@np.arange(n)).astype(int)
print(1000*n-value(m.objective)) # 最適化後の総移動時間
>>>
Wall time: 488 ms
1973.67977958
2/3以下になりました。
分布を出してみましょう。
python
plt.title('通勤時間分布')
plt.hist([1000-g.edge[i][s[i]+100]['weight'] for i in range(n)], alpha=0.5, range=(0,180), bins=60)
plt.hist([1000-g.edge[i][i+100]['weight'] for i in range(n)], alpha=0.5, range=(0,180), bins=60);
ちなみに、NetworkXでも解けますが、遅いです。
python
%time r = nx.max_weight_matching(g)
print(1000*n-sum(g.edge[i][r[i]]['weight'] for i in range(n)))
>>>
Wall time: 2.73 s
1973.67977958
以上