이게 뭐야

4월부터 새롭게 하숙하고, 통근·통학을 시작한 분도 많을지도 모릅니다.
도쿄라면 좀처럼 근처에 사는 것도 어렵습니다.

그런데, 모두가 사는 곳을 셔플 해, 자유롭게 이사 가능하게 한 경우, 총 통근 시간을 최소화하면 어떻게 될까요?
(근무지도 주택의 위치도 같고, 사는 사람을 바꾸는 것으로 합니다)

풀어봐

통근 시간의 분포입니다만, 통근에 관한 설문 조사 결과 의 전철 통근을 보면, 평균 33.23분, 표준 편차 21.80분으로 있으므로, 그것을 사용합니다. 분포는 적당히 로그 정규 분포로 합시다.

주택은 반경 60분의 원 안에 균일하게 무작위로 만듭니다. 근무지는, 거기로부터 방향은 랜덤으로, 상기의 거리 떨어진 곳으로 합니다.

파이썬으로 보자.

파이썬

%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

파이썬

plt.title('住宅')
plt.scatter(xs,ys);

파이썬

plt.title('勤務地')
plt.scatter(xt,yt);

파이썬

print(dist2.sum()) # 元の総移動時間
>>>
3323.3549788106984

최적화해 봅시다.

파이썬

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 이하가 되었습니다.
분포를 내보자.

파이썬

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에서도 풀 수 있지만, 느립니다.

파이썬

%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

이상