p-center 問題入門¶
1. はじめに¶
多くの自治体では、消防署や避難所といった拠点をなるべく「どの住民からも遠過ぎない」場所に建てようとする。そこで役に立つのがp-center問題である。このモデルでは、
- 対象地域上の住民が集まる点（需要点）から最寄りの施設まで徒歩や車で行く距離を求める（需要点は複数存在する）。
- その距離の中で「いちばん長いもの」をサービス距離と呼ぶ。
- このサービス距離をできるだけ小さくするように、
p箇所だけ施設を置く場所を選ぶ。
つまり、いちばん困っている人をいちばん助けるという考え方で場所を決めるわけである。
2. 前提知識と学習目標¶
本資料を読む前に知っておくと良いことは次のとおりである。
- 数学(高校 1–2 年程度): 比や比例、最大値、最小値、場合分けで考える方法。
- プログラミング:
if文で条件分岐し、for文で同じ処理を繰り返し、listでデータをまとめて扱える。 - Google Colab: Python コードセルを実行できる。
学習目標¶
- p-center 問題が「最悪の距離」を小さくするモデルであることを説明できる。
- 具体例から数式モデルを組み立てられる。
- PuLP でモデルを作り、答えを得て、図で示せる。
3. 問題設定と定式化¶
3.1 問題の具体例¶
消防署を例に具体的な数字を使って説明する。
- 候補点: 消防署を建設できそうな場所 4 か所 $I=\{1,\dots,4\}$。
- 需要点: 町や丁などの中心地点 8 か所 $J=\{1,\dots,8\}$。
- 設置できる消防署の数: $p=2$。
- 各候補点と各需要点との距離は、直線距離とする。
3.2 数学的な定式化¶
3.2.1 定数および決定変数の定義¶
| 名前 | 説明 |
|---|---|
| $d_{ij}$ | 候補点 $i$ から需要点 $j$ までの距離 |
| $p$ | 設置する施設の数（例では 2） |
| $x_{ij}$ | 需要点 $j$ を候補点 $i$ へ割り当てるなら 1、そうでなければ 0 |
| $y_i$ | 候補点 $i$ に施設を建てるなら 1、そうでなければ 0 |
| $R$ | すべての需要点が最寄り施設に行くときの最長距離 |
3.2.2 数理モデル¶
目的は $R$ を最小にすることである。
$$ \min R $$
制約は次の 4 種である。
- 必ずどこかへ割り当てる
$$\sum_{i \in I} x_{ij} = 1 \quad (\forall j \in J)$$ - 割当は建てた場所にだけ
$$x_{ij} \le y_i \quad (\forall i \in I,\; j \in J)$$ - 建設数は $p$ 個だけ
$$\sum_{i \in I} y_i = p$$ - 最長距離の定義
$$d_{ij}\,x_{ij} \le R \quad (\forall i \in I,\; j \in J)$$
$x_{ij}, y_i$ は 0 か 1 だけを取る。$R$ は 0 以上の実数である。
3.2.3 目的関数の意味¶
$$\min R$$
この式は、需要点それぞれが「いちばん近い施設」に行く距離を計算したあと、その中でいちばん大きいものを $R$ とし、それを小さくしようとするものである。例えば、最寄り施設まで 400 m かかる家が 1 軒だけあり、他は 300 m 以内なら $R=400$ である。$R$ が小さいほど、遠い人でも短く移動できるので、公平性を守る配置になる。施設を増やせば $R$ は小さくしやすくなるが、$p$ が固定なら、限られた数で最善を尽くすのが p-center 問題である。
3.2.4 制約条件の意味¶
必ずどこかへ割り当てる
$$\sum_{i \in I} x_{ij} = 1$$ 需要点 $j$ は必ず 1 か所の施設に割り当てられる。0 や 2 になることはない。これにより「どの需要点に住む住民でもサービスを受けられる」ことが保証される。割当は建てた場所にだけ
$$x_{ij} \le y_i$$ 候補点 $i$ に施設を建てたとき（$y_i=1$）だけ、需要点 $j$ をそこに割り当てる（$x_{ij}=1$）ことが許される。$y_i=0$ のときは $x_{ij}$ も 0 しか選べない。建設数は $p$ 個だけ
$$\sum_{i \in I} y_i = p$$ 予算や人員の都合で、建てられる施設は合計で $p$ 箇所ちょうどである。多過ぎても少な過ぎてもダメ、というルールを一行で書いている。最長距離の定義
$$d_{ij}\,x_{ij} \le R$$ 需要点 $j$ が候補点 $i$ に割り当てられた場合だけ（$x_{ij}=1$）、距離 $d_{ij}$ が $R$ 以下でなければならない。つまり、すべての割り当てについて $R$ が上限になるように選ばれ、結果として $R$ は「最長距離」を表す変数として働く。
4. LPファイルとの対応関係¶
4.1 LPファイルとは¶
LP ファイルは、最適化ソルバー（Gurobi や CBC など）が読めるテキスト形式である。モデル名、目的関数、制約、変数の種類を順番に並べたものに過ぎない。PuLP では writeLP() で簡単に出力できる。
4.2 p-center モデルと LP ファイルの対応¶
例として、需要点 3 か所、候補点 2 か所、$p=1$ の小規模な問題を出力すると、次のようになる。
lp
Minimize
obj: R
Subject To
c_assign_1: x_1_1 + x_2_1 = 1
c_assign_2: x_1_2 + x_2_2 = 1
c_assign_3: x_1_3 + x_2_3 = 1
c_link_1_1: x_1_1 - y_1 <= 0
...
c_p: y_1 + y_2 = 1
c_R_1_1: 250 x_1_1 - R <= 0
...
Binary
x_1_1 x_1_2 ... y_1 y_2
End
ここで c_assign_1 などは制約の名前、obj が目的関数である。式自体は 3.2 節と同じ意味をもつ。
5. Python + PuLP による実装¶
5.1 準備¶
以下にp-center 問題の実装例を示す。 問題は候補点4 点、需要点8 点 とかなり小規模な問題である。 頑張れば手でも解けるので、ぜひ試してみてほしい。
なおGoogle Colab で実行する際には、必ず最初に以下のプログラム2行目（# !pip install pulp matplotlib japanize-matplotlib）のコメントを外して実行すること。
In [1]:
# Google Colab の場合は↓のコメントを外して実行: ライブラリインストール（1 回だけ）
# !pip install pulp matplotlib japanize-matplotlib
In [2]:
# ライブラリ読み込み
import pulp # モデル作成
import matplotlib.pyplot as plt # 可視化
# ダミーデータ作成（候補点 4 点、需要点 8 点）
coords_fac = [(0,0), (10,0), (0,10), (10,10)] # 候補点座標
coords_dem = [(2,1), (3,7), (6,2), (7,8), (1,9), (9,1), (5,5), (8,4)] # 需要点座標
I = range(len(coords_fac)) # 候補点集合
J = range(len(coords_dem)) # 需要点集合
p = 2 # 配置できる施設の数（この例では 2 か所）
print(f"候補点数: {len(I)}, 需要点数: {len(J)}")
候補点数: 4, 需要点数: 8
In [3]:
import math
import numpy as np
dist = np.zeros((len(I), len(J))) # 距離行列（例：0行3列なら候補点0 から需要点3 までの距離）
for i in I:
for j in J:
# 距離計算（ユークリッド距離）
dist[i][j] = math.sqrt((coords_fac[i][0] - coords_dem[j][0])**2 + (coords_fac[i][1] - coords_dem[j][1])**2)
### 距離行列表示ここから
# 以下は候補点から需要点までの距離するためのプログラム
# わからなければ読み飛ばしてOK
import unicodedata
w = lambda t: sum(2 if unicodedata.east_asian_width(c) in 'FW' else 1 for c in t)
width = w("候補点0 から需要点j までの距離")
print(f"{' '*(width-11)}需要点j =", end=" ")
for j in J:
print(j, end=" " * (len(f"{dist[0][j]:.2f} km,")))
print()
for i in I:
print(f"候補点{i} から需要点j までの距離", end=" ")
for j in J:
print(f"{dist[i][j]:.2f} km", end=", ")
print()
### 距離行列表示ここまで
### 地図表示ここから
legend_flg = True
for i, (x_f, y_f) in enumerate(coords_fac):
if legend_flg:
plt.scatter(x_f, y_f, c='gray', s=100, label=f"selected facilities")
legend_flg = False
else:
plt.scatter(x_f, y_f, c='gray', s=100)
plt.text(x_f+0.3, y_f, f"cand{i}")
legend_flg = True
for j, (x_d, y_d) in enumerate(coords_dem):
if legend_flg:
plt.scatter(x_d, y_d, c='blue', s=40, label=f"demand points")
legend_flg = False
else:
plt.scatter(x_d, y_d, c='blue', s=40)
plt.text(x_d+0.3, y_d, f"dem{j}")
需要点j = 0 1 2 3 4 5 6 7 候補点0 から需要点j までの距離 2.24 km, 7.62 km, 6.32 km, 10.63 km, 9.06 km, 9.06 km, 7.07 km, 8.94 km, 候補点1 から需要点j までの距離 8.06 km, 9.90 km, 4.47 km, 8.54 km, 12.73 km, 1.41 km, 7.07 km, 4.47 km, 候補点2 から需要点j までの距離 9.22 km, 4.24 km, 10.00 km, 7.28 km, 1.41 km, 12.73 km, 7.07 km, 10.00 km, 候補点3 から需要点j までの距離 12.04 km, 7.62 km, 8.94 km, 3.61 km, 9.06 km, 9.06 km, 7.07 km, 6.32 km,
灰点（candX）は候補点、青点（demX）は需要点を表しており、Xは番号である。
この問題では図の四つ角にある候補点cand0~cand3 の中から、この地域住民（dem0~dem7 に住んでる人たち）にとって最も公平な点を二つ選んで施設を配置する。
なお、「公平」とはモデルの説明で述べた通り、「施設から最も遠い需要点」が「最も短くなるような候補点」を選ぶことを指す。
5.2 モデル構築¶
In [4]:
# モデル宣言
prob = pulp.LpProblem("p_center", pulp.LpMinimize)
# 変数
x = pulp.LpVariable.dicts('x', (I, J), 0, 1, cat='Binary')
y = pulp.LpVariable.dicts('y', I, 0, 1, cat='Binary')
R = pulp.LpVariable('R', lowBound=0)
# 目的関数
prob += R
# 制約 1: 割当必須
for j in J:
prob += pulp.lpSum(x[i][j] for i in I) == 1
# 制約 2: 割当は開設先へ
for i in I:
for j in J:
prob += x[i][j] <= y[i]
# 制約 3: 設置数
prob += pulp.lpSum(y[i] for i in I) == p
# 制約 4: R の定義（ユークリッド距離）
for i in I:
for j in J:
prob += dist[i][j] * x[i][j] <= R
print(prob)
p_center: MINIMIZE 1*R + 0 SUBJECT TO _C1: x_0_0 + x_1_0 + x_2_0 + x_3_0 = 1 _C2: x_0_1 + x_1_1 + x_2_1 + x_3_1 = 1 _C3: x_0_2 + x_1_2 + x_2_2 + x_3_2 = 1 _C4: x_0_3 + x_1_3 + x_2_3 + x_3_3 = 1 _C5: x_0_4 + x_1_4 + x_2_4 + x_3_4 = 1 _C6: x_0_5 + x_1_5 + x_2_5 + x_3_5 = 1 _C7: x_0_6 + x_1_6 + x_2_6 + x_3_6 = 1 _C8: x_0_7 + x_1_7 + x_2_7 + x_3_7 = 1 _C9: x_0_0 - y_0 <= 0 _C10: x_0_1 - y_0 <= 0 _C11: x_0_2 - y_0 <= 0 _C12: x_0_3 - y_0 <= 0 _C13: x_0_4 - y_0 <= 0 _C14: x_0_5 - y_0 <= 0 _C15: x_0_6 - y_0 <= 0 _C16: x_0_7 - y_0 <= 0 _C17: x_1_0 - y_1 <= 0 _C18: x_1_1 - y_1 <= 0 _C19: x_1_2 - y_1 <= 0 _C20: x_1_3 - y_1 <= 0 _C21: x_1_4 - y_1 <= 0 _C22: x_1_5 - y_1 <= 0 _C23: x_1_6 - y_1 <= 0 _C24: x_1_7 - y_1 <= 0 _C25: x_2_0 - y_2 <= 0 _C26: x_2_1 - y_2 <= 0 _C27: x_2_2 - y_2 <= 0 _C28: x_2_3 - y_2 <= 0 _C29: x_2_4 - y_2 <= 0 _C30: x_2_5 - y_2 <= 0 _C31: x_2_6 - y_2 <= 0 _C32: x_2_7 - y_2 <= 0 _C33: x_3_0 - y_3 <= 0 _C34: x_3_1 - y_3 <= 0 _C35: x_3_2 - y_3 <= 0 _C36: x_3_3 - y_3 <= 0 _C37: x_3_4 - y_3 <= 0 _C38: x_3_5 - y_3 <= 0 _C39: x_3_6 - y_3 <= 0 _C40: x_3_7 - y_3 <= 0 _C41: y_0 + y_1 + y_2 + y_3 = 2 _C42: - R + 2.2360679775 x_0_0 <= 0 _C43: - R + 7.61577310586 x_0_1 <= 0 _C44: - R + 6.32455532034 x_0_2 <= 0 _C45: - R + 10.6301458127 x_0_3 <= 0 _C46: - R + 9.05538513814 x_0_4 <= 0 _C47: - R + 9.05538513814 x_0_5 <= 0 _C48: - R + 7.07106781187 x_0_6 <= 0 _C49: - R + 8.94427191 x_0_7 <= 0 _C50: - R + 8.0622577483 x_1_0 <= 0 _C51: - R + 9.89949493661 x_1_1 <= 0 _C52: - R + 4.472135955 x_1_2 <= 0 _C53: - R + 8.54400374532 x_1_3 <= 0 _C54: - R + 12.7279220614 x_1_4 <= 0 _C55: - R + 1.41421356237 x_1_5 <= 0 _C56: - R + 7.07106781187 x_1_6 <= 0 _C57: - R + 4.472135955 x_1_7 <= 0 _C58: - R + 9.21954445729 x_2_0 <= 0 _C59: - R + 4.24264068712 x_2_1 <= 0 _C60: - R + 10 x_2_2 <= 0 _C61: - R + 7.28010988928 x_2_3 <= 0 _C62: - R + 1.41421356237 x_2_4 <= 0 _C63: - R + 12.7279220614 x_2_5 <= 0 _C64: - R + 7.07106781187 x_2_6 <= 0 _C65: - R + 10 x_2_7 <= 0 _C66: - R + 12.0415945788 x_3_0 <= 0 _C67: - R + 7.61577310586 x_3_1 <= 0 _C68: - R + 8.94427191 x_3_2 <= 0 _C69: - R + 3.60555127546 x_3_3 <= 0 _C70: - R + 9.05538513814 x_3_4 <= 0 _C71: - R + 9.05538513814 x_3_5 <= 0 _C72: - R + 7.07106781187 x_3_6 <= 0 _C73: - R + 6.32455532034 x_3_7 <= 0 VARIABLES R Continuous 0 <= x_0_0 <= 1 Integer 0 <= x_0_1 <= 1 Integer 0 <= x_0_2 <= 1 Integer 0 <= x_0_3 <= 1 Integer 0 <= x_0_4 <= 1 Integer 0 <= x_0_5 <= 1 Integer 0 <= x_0_6 <= 1 Integer 0 <= x_0_7 <= 1 Integer 0 <= x_1_0 <= 1 Integer 0 <= x_1_1 <= 1 Integer 0 <= x_1_2 <= 1 Integer 0 <= x_1_3 <= 1 Integer 0 <= x_1_4 <= 1 Integer 0 <= x_1_5 <= 1 Integer 0 <= x_1_6 <= 1 Integer 0 <= x_1_7 <= 1 Integer 0 <= x_2_0 <= 1 Integer 0 <= x_2_1 <= 1 Integer 0 <= x_2_2 <= 1 Integer 0 <= x_2_3 <= 1 Integer 0 <= x_2_4 <= 1 Integer 0 <= x_2_5 <= 1 Integer 0 <= x_2_6 <= 1 Integer 0 <= x_2_7 <= 1 Integer 0 <= x_3_0 <= 1 Integer 0 <= x_3_1 <= 1 Integer 0 <= x_3_2 <= 1 Integer 0 <= x_3_3 <= 1 Integer 0 <= x_3_4 <= 1 Integer 0 <= x_3_5 <= 1 Integer 0 <= x_3_6 <= 1 Integer 0 <= x_3_7 <= 1 Integer 0 <= y_0 <= 1 Integer 0 <= y_1 <= 1 Integer 0 <= y_2 <= 1 Integer 0 <= y_3 <= 1 Integer
表示されているのはpulp によって生成されたlp ファイルの中身である。
5.3 求解および結果の解釈¶
In [5]:
# CBC （scip のようなもの）ソルバーで解く
prob.solve(pulp.PULP_CBC_CMD(msg=1))
print(f"最長距離 R* = {pulp.value(R):.2f} km")
# 開設した施設
open_sites = [i for i in I if pulp.value(y[i]) > 0.5]
print("開設候補点番号:", open_sites)
Welcome to the CBC MILP Solver Version: 2.10.10 Build Date: Sep 26 2023 command line - cbc /tmp/9c8588f56b8e42658fca6b8f9be54300-pulp.mps -timeMode elapsed -branch -printingOptions all -solution /tmp/9c8588f56b8e42658fca6b8f9be54300-pulp.sol (default strategy 1) At line 2 NAME MODEL At line 3 ROWS At line 78 COLUMNS At line 316 RHS At line 390 BOUNDS At line 427 ENDATA Problem MODEL has 73 rows, 37 columns and 164 elements Coin0008I MODEL read with 0 errors Option for timeMode changed from cpu to elapsed Continuous objective value is 1.76777 - 0.00 seconds Cgl0004I processed model has 73 rows, 37 columns (36 integer (36 of which binary)) and 164 elements Cbc0038I Initial state - 35 integers unsatisfied sum - 9.41862 Cbc0038I Pass 1: suminf. 0.00000 (0) obj. 12.7279 iterations 38 Cbc0038I Solution found of 12.7279 Cbc0038I Relaxing continuous gives 12.7279 Cbc0038I Before mini branch and bound, 1 integers at bound fixed and 0 continuous Cbc0038I Full problem 73 rows 37 columns, reduced to 71 rows 36 columns Cbc0038I Mini branch and bound improved solution from 12.7279 to 8.06226 (0.01 seconds) Cbc0038I Round again with cutoff of 7.4328 Cbc0038I Pass 2: suminf. 1.47819 (8) obj. 7.4328 iterations 4 Cbc0038I Pass 3: suminf. 1.47819 (8) obj. 7.4328 iterations 0 Cbc0038I Pass 4: suminf. 1.45828 (6) obj. 7.4328 iterations 12 Cbc0038I Pass 5: suminf. 1.14021 (6) obj. 7.4328 iterations 5 Cbc0038I Pass 6: suminf. 2.39273 (15) obj. 7.4328 iterations 11 Cbc0038I Pass 7: suminf. 1.14021 (6) obj. 7.4328 iterations 9 Cbc0038I Pass 8: suminf. 1.45828 (6) obj. 7.4328 iterations 9 Cbc0038I Pass 9: suminf. 1.45828 (6) obj. 7.4328 iterations 3 Cbc0038I Pass 10: suminf. 1.45828 (6) obj. 7.4328 iterations 3 Cbc0038I Pass 11: suminf. 3.59165 (15) obj. 7.4328 iterations 11 Cbc0038I Pass 12: suminf. 1.61443 (8) obj. 7.4328 iterations 13 Cbc0038I Pass 13: suminf. 1.61443 (8) obj. 7.4328 iterations 6 Cbc0038I Pass 14: suminf. 1.61443 (8) obj. 7.4328 iterations 4 Cbc0038I Pass 15: suminf. 1.14021 (6) obj. 7.4328 iterations 6 Cbc0038I Pass 16: suminf. 1.45828 (6) obj. 7.4328 iterations 9 Cbc0038I Pass 17: suminf. 4.52678 (23) obj. 7.4328 iterations 16 Cbc0038I Pass 18: suminf. 1.79625 (8) obj. 7.4328 iterations 9 Cbc0038I Pass 19: suminf. 1.14021 (6) obj. 7.4328 iterations 7 Cbc0038I Pass 20: suminf. 1.45828 (6) obj. 7.4328 iterations 9 Cbc0038I Pass 21: suminf. 1.45828 (6) obj. 7.4328 iterations 8 Cbc0038I Pass 22: suminf. 2.03224 (17) obj. 7.4328 iterations 15 Cbc0038I Pass 23: suminf. 1.84003 (8) obj. 7.4328 iterations 12 Cbc0038I Pass 24: suminf. 1.11683 (6) obj. 7.4328 iterations 10 Cbc0038I Pass 25: suminf. 1.48166 (6) obj. 7.4328 iterations 6 Cbc0038I Pass 26: suminf. 3.64472 (15) obj. 7.4328 iterations 12 Cbc0038I Pass 27: suminf. 1.84003 (8) obj. 7.4328 iterations 9 Cbc0038I Pass 28: suminf. 1.11683 (6) obj. 7.4328 iterations 10 Cbc0038I Pass 29: suminf. 1.48166 (6) obj. 7.4328 iterations 6 Cbc0038I Pass 30: suminf. 1.63781 (8) obj. 7.4328 iterations 2 Cbc0038I Pass 31: suminf. 1.63781 (8) obj. 7.4328 iterations 0 Cbc0038I No solution found this major pass Cbc0038I Before mini branch and bound, 0 integers at bound fixed and 0 continuous Cbc0038I Full problem 73 rows 37 columns, reduced to 73 rows 37 columns Cbc0038I Mini branch and bound did not improve solution (0.01 seconds) Cbc0038I After 0.01 seconds - Feasibility pump exiting with objective of 8.06226 - took 0.00 seconds Cbc0012I Integer solution of 8.0622577 found by feasibility pump after 0 iterations and 0 nodes (0.01 seconds) Cbc0006I The LP relaxation is infeasible or too expensive Cbc0013I At root node, 0 cuts changed objective from 1.767767 to 1.767767 in 1 passes Cbc0014I Cut generator 0 (Probing) - 1 row cuts average 0.0 elements, 3 column cuts (3 active) in 0.000 seconds - new frequency is 1 Cbc0014I Cut generator 1 (Gomory) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 2 (Knapsack) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 3 (Clique) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 4 (MixedIntegerRounding2) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 5 (FlowCover) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 6 (TwoMirCuts) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0014I Cut generator 7 (ZeroHalf) - 0 row cuts average 0.0 elements, 0 column cuts (0 active) in 0.000 seconds - new frequency is -100 Cbc0001I Search completed - best objective 8.062257748299, took 0 iterations and 0 nodes (0.01 seconds) Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost Cuts at root node changed objective from 1.76777 to 1.76777 Probing was tried 1 times and created 4 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds) Result - Optimal solution found Objective value: 8.06225775 Enumerated nodes: 0 Total iterations: 0 Time (CPU seconds): 0.01 Time (Wallclock seconds): 0.01 Option for printingOptions changed from normal to all Total time (CPU seconds): 0.01 (Wallclock seconds): 0.01
最長距離 R* = 8.06 km 開設候補点番号: [1, 2]
ここでは5.2 で定式化したp-center 問題をコンピュータに解かせている。
Result - Optimal solution found となっているので最適解が見つかったようである。
最後の2行のうち、1行目「最長距離 R*」は「配置された施設から最も離れた需要点までの距離（最適値）が8.06 km だった」という意味である。
2行目は「開設候補点番号」は問題を解いて施設が配置された候補点（最適解）が「候補点1, 2」であることを意味している。
つまり、四つの候補点のうち最適解として選ばれた1, 2 から最も遠い需要点までの距離は高々8.06 km でこの8.06 km という距離は最悪ケースの中で最良であることがわかる
（適当に配置点を選ぶと最も遠い需要点までの距離は必ず8.06 km 以上になってしまう）。
5.4 可視化¶
In [6]:
import japanize_matplotlib
# 座標を散布図で表示
assign_mat = np.array([[pulp.value(x[i][j]) for i in I] for j in J])
assigned = np.where(assign_mat)
legend_flg = True
for asgn_i, asgn_j in zip(assigned[1], assigned[0]):
x_f, y_f = coords_fac[asgn_i]
x_d, y_d = coords_dem[asgn_j]
if legend_flg:
plt.plot([x_d, x_f], [y_d, y_f], c='black', lw=1, alpha=0.6, label="assignment")
legend_flg = False
else:
plt.plot([x_d, x_f], [y_d, y_f], c='black', lw=1, alpha=0.6)
legend_flg = True
for i, (x_f, y_f) in enumerate(coords_fac):
if i in open_sites:
if legend_flg:
plt.scatter(x_f, y_f, c='red', s=100, label=f"selected facilities")
legend_flg = False
else:
plt.scatter(x_f, y_f, c='red', s=100)
else:
plt.scatter(x_f, y_f, c='gray', s=100)
plt.text(x_f+0.3, y_f, f"cand{i}")
legend_flg = True
for j, (x_d, y_d) in enumerate(coords_dem):
if legend_flg:
plt.scatter(x_d, y_d, c='blue', s=40, label=f"demand points")
legend_flg = False
else:
plt.scatter(x_d, y_d, c='blue', s=40)
plt.text(x_d+0.3, y_d, f"dem{j}")
plt.title('p-center 配置結果')
plt.legend()
plt.show()
赤点は施設の建設地点として選ばれた候補点を意味する。 黒線は実際に配置された施設がどの需要点に割り当てられているかを示している。 上の例だと、施設は候補点1 と2 に配置され、候補点1 の施設は需要点0, 2, 5, 6, 7 に割り当てられ、候補点2 の施設は需要点1, 3, 4 に割り当てられていることがわかる。
6. 注意事項¶
6.1 プログラムに関して¶
- このプログラムでは候補点と需要点との距離を、直線距離で計算している。実際の問題では必要に応じてGoogle map で算出した移動時間などに置き換えること。
- 候補点や需要点、施設配置数$p$ を増やすと計算時間が急に増える場合がある。計算不能になったら指導教員に相談すること。
- コードを改変するときは、まず小さなデータで動作を確認してから大きな地図データに適用するとエラー箇所を特定しやすくなる。
- 今回の小規模問題でp-center 問題の目的や意味が理解できてから実際（例えば習志野市）の問題に取り組むこと。
- 現実のデータは国土交通省や統計局が提供している地区データや人口データを利用することで取得することができる。
- 習志野市のデータが欲しいゼミ生は別途指導教員に相談すること。
6.2 自分で問題を作る場合には¶
5.1 準備 に記載のプログラムのうち、以下の箇所を変更すると自分でダミー問題を作ることができる。
# ダミーデータ作成（候補点 4 点、需要点 8 点）
coords_fac = [(0,0), (10,0), (0,10), (10,10)] # 候補点座標
coords_dem = [(2,1), (3,7), (6,2), (7,8), (1,9), (9,1), (5,5), (8,4)] # 需要点座標
coords_fac に格納されている(0, 0) は(x座標、y座標) という意味である。
coords_fac の(x座標、y座標) を増やせば候補点を増やすことができる。
coords_dem の(x座標、y座標) を増やせば需要点を増やすことができる。
また、(x座標、y座標) の座標をGoogle map などで調べた地図上の(経度, 緯度) に置き換えれば、現実の座標で問題を作ることもできる。
なお、一般的には(緯度, 経度) という順番で表記するが、実装上は経度をx 座標、緯度をy 座標としたほうが可視化の際に混乱しない。
色々なパターンで点を追加して試してみるとよい。
また、候補点および需要点を増やしたら、施設配置数p も変更して問題を作ってみると良いだろう。
6.3 定式化が理解できないときは¶
数理最適化問題の定式化に関して、慣れないうちは理解が難しいと感じるかもしれない。 なかなか理解が進まない人の特徴として数式だけで考えようとすることが挙げられる。 数式だけで考えるのではなく、実際の問題をイメージしながら定式化を読み解くことが重要である。 そのため、定式化されたモデルを理解するには、小さな問題を手作業で解いてみることが大切である。 p-center 問題であれば、例えば2つの候補点、4つの需要点を持つ問題を考えてみて、ノートに実際に点を書き出し、距離行列を作ってみること。 その上で、用意した小規模問題のデータを使って、定式化された数式を展開し、実際に最適解と最適値を求めてみること。 このとき大事なのは、プログラムは一切使用せず、全て手作業で行うことである。 （問題が小さければ手作業で解いた問題が本当に最適であるか、制約条件に値を代入することで簡単に確認することができる）
小規模問題の作成と手作業による最適化が終わったら、今度はその小規模問題をプログラムで解いてみること。 このとき、手作業で解いた結果とプログラムの結果が一致することを確認すること。 一致していれば、おそらく大きな問題に対しても正しく定式化できているし、一致していなければプログラムが間違っている。 プログラム上のバグを見つけるためには、print文を使って、細かく各変数（ここでいう変数とはプログラム上の変数のことである）の値をチェックし、時間をかけることが大切である。
これらの手続きを踏んでも理解できない場合には、指導教員まで相談すること。
7. 卒業研究概要への手引き¶
- このテーマを卒業研究のテーマにしたい場合は、ここの対応をすること(必ず事前に相談すること)。
- 以下、すべて含めて、A4用紙2枚にまとめること。
- 厳密に守る必要はないが、文章量の比率を目安にすること(15%とは、A4用紙2枚分の15%という意味)。
1章 はじめに(文章量の比率: 15%)¶
- 施設配置問題とは何か、p-center問題以外にどのような施設配置問題があるのかまとめること。
- p-center問題に関する先行研究を紹介すること。
- 施設配置問題を解く際に使われる整数計画法とは何か、まとめること。
2章 p-center問題(文章量の比率: 30%)¶
- p-center問題の数理モデルについて、説明すること。最も遠い需要点までの距離を最小化するという、p-median問題との目的関数の違いも説明すること。
3章 実験(文章量の比率: 40%)¶
3.1節 概要¶
- 本資料のダミーデータをベースに、実際に新たな地域(需要点と候補点)を作成して問題を解き、図示すること。
- 施設数 $p$ を変えて解の変化を確認すること。
3.2節 結果と考察¶
- サービス距離(最も遠い需要点までの距離)がいくつになったか示すこと。
- そのサービス距離が何を意味しているか説明すること。
- 図示した画像から、どのようなことが言えるか、考え説明すること。
4章 おわりに(文章量の比率: 15%)¶
- 今回の実験に対する感想を記載すること。
- 例えば、直感で施設を配置するのと数理モデルにより配置するのでは、どちらがよさそうか、また良い理由を記載すること。
参考文献¶
参考にした資料を、2〜3件記載すること。以下、書き方の例である。
- [1] 柿本ほか, XXXに関する分析, XXX学会論文誌, 2020.
- [2] XXXに関する情報, http://xxx.ddd.ttt.com, 2020年4月20日閲覧
本文中で引用する場合は「柿本らはXXXを実施している [1]。また、〜」のように、どこで引用したのか明白にすること。