PyLearn MAP 04: 遺伝的アルゴリズムによる最適施設配置の近似解¶
日本大学生産工学部 マネジメント工学科
知能情報システム工学研究室
大前佑斗
PyLearnMAP03では、現在、習志野市にある5箇所の消防署がどこに配置されているのかプロットする方法、市民1人当たりとどれくらい距離が離れているのか（利便性）、算出する方法と学びました。それで、現在の配置は1人当たりの平均として、828m離れていることがわかりました。ここで、もっと良い配置はないのかと疑問が湧きます。
このような問題、つまり、施設の最適配置は、オペレーションズリサーチ（OR）の研究領域になります。その中でも特に、Facility Location Problem (FLP) と呼びます。ここでは、最も良い施設配置に近い場所、すなわち近似解を求める手法、遺伝的アルゴリズムを活用して、習志野市におけるより良い施設配置を探す方法について解説します。
0. 前回の復習¶
まず、前回の復習として、現在の習志野市の消防署5箇所と、500mメッシュ地点の人口をプロットします。前回のコードそのまんまですので、特に解説はしません。
import folium
import json
# 市町村の境界データの読み込み
f = open("N03-19_12_190101.geojson", 'r', encoding='utf-8_sig')
geojson = json.load(f)
f.close()
# 中心取得
import geocoder as ge
cent = ge.arcgis("習志野市")
cent = cent.latlng
# 習志野市が格納されている区画のIDを取得
Target = "習志野市"
for i in range(len(geojson["features"])):
txt = geojson["features"][i]["properties"]["N03_004"]
if txt == Target:
GetID = i
break
# 習志野市の人口データ取得
import shapefile
src = shapefile.Reader('Mesh4_POP_12.shx',encoding='SHIFT-JIS')
TargetID = 12216 # 習志野市の市町村コード
Narashino = []
for i in range(0, len(src)):
if src.record(i)[1] == TargetID:
Narashino.append(i)
Place = []
MeshCode = []
Population = []
# (なぜか0番目に埼玉が混ざっているので1から)
for i in range(1, len(src)): # len(src)は千葉県内のメッシュの数
# 習志野市であれば、経度緯度、人口をappend
for j in range(0, len(Narashino)):
if Narashino[j] == i:
# 座標の取得
a = src.shape(i).points[0][0]
b = src.shape(i).points[0][1]
Place.append([float(b), float(a)]) # 経度と緯度を逆にし、Foliumの順に
# メッシュコードの取得
MeshCode.append(str(src.record(i)[0]))
# 2010年の総人口を取得
Population.append(src.record(i)[2])
# 地図の描画
# 500mメッシュ地点のマーカを設定
from folium.features import DivIcon
# 赤いボックス
def MakeBox_red(in_text):
out_text = '''
<div style="text-align: center;
font-size: 5pt;
color: white;
padding-top: 0px;
height: 10px;
width: 23px;
background: red;
border: 0.5px solid;
border-radius: 30%;">'''
out_text = out_text+ in_text + '</div>'
return out_text
# 青いボックス
def MakeBox_blue(in_text):
out_text = '''
<div style="text-align: center;
font-size: 5pt;
color: white;
padding-top: 0px;
height: 10px;
width: 23px;
background: blue;
border: 0.5px solid;
border-radius: 30%;">'''
out_text = out_text+ in_text + '</div>'
return out_text
FireName = [
"千葉県習志野市鷺沼２丁目１−４３",
"千葉県習志野市東習志野２丁目２−１５",
"千葉県習志野市藤崎６丁目２０−１１",
"千葉県習志野市奏の杜２丁目１３−１",
"千葉県習志野市秋津３丁目７−１"]
FirePoint = []
for i in range(len(FireName)):
FirePoint.append(ge.arcgis(FireName[i]).latlng)
from geopy import distance
import numpy as np
# FirePoint: リスト化された、消防署の経度・緯度
# Place: リスト化された、メッシュ地点の経度・緯度
def DecisionPoint(FirePoint, Place):
# 以下の処理を全メッシュ地点に行う
MinPoint = []
for i in range(len(Place)):
# メッシュ1地点と、消防署全地点の距離を求める
temp_list = []
for j in range(len(FirePoint)):
temp_list.append(distance.distance(Place[i], FirePoint[j]).km)
# 距離が最小の消防署のIDを取得
MinPoint.append(np.argmin(temp_list))
return MinPoint
# メッシュ地点に対する、利用する消防署IDを取得
MinPoint = DecisionPoint(FirePoint, Place)
# 習志野市を中心とするマップの作成
m6 = folium.Map(location=cent, zoom_start=13)
# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
fill_color='YlGn', # 塗りつぶしの色
fill_opacity=0.1, # 塗りつぶしの濃さ
line_opacity=0.9, # 境界の色の濃さ
).add_to(m6)
# 500mメッシュ地点のマーカを設定
MeanPop = np.mean(Population)
for i in range(0, len(Place)):
if Population[i] > MeanPop: # 1万人オーバーなら、青
html_text = str(int(Population[i]))
folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m6)
else:
html_text = str(int(Population[i]))
folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m6)
# 各消防署のカバー人口を算出
CoverPop = []
for i in range(len(FirePoint)):
s = 0
for j in range(len(Population)):
if MinPoint[j] == i: # そのメッシュ地点に対応している消防署のみをカウント
s = s + Population[j]
CoverPop.append(s) # 消防署iが対応する総人数をappend
# 消防署をプロット
for i in range(len(FireName)):
folium.Marker(location = FirePoint[i],
popup = "<p>Fire station: "+ str(i)+"</p><p>Cover: " + str(int(CoverPop[i])) + " persons</p>").add_to(m6)
# 消防署と500mメッシュのラインを表示
for i in range(len(Place)):
folium.PolyLine(locations=[FirePoint[MinPoint[i]], Place[i]], weight=2, color='black').add_to(m6)
# 凡例はhtmlで入れるしかない模様
legend_html = '''
<div style="position: fixed;
top: 20px; left: 50px; width: 210px; height: 67px; Background: white;
border:2px solid grey; z-index:9999; font-size:14px; padding-left: 5px">
Blue box: High population area.<br>
Red box: Low population area.<br>
Marker: Fire station
</div>
'''
legend_html2 = '''
<div style="position: fixed;
top: 1px; left: 50px; width: 260px; height: 67px; Background: transparent;;
border:0px solid grey; z-index:9999; font-size:14px; padding-left: 5px">
Made by Yuto Omae at Nihon Univ.
</div>
'''
m6.get_root().html.add_child(folium.Element(legend_html))
m6.get_root().html.add_child(folium.Element(legend_html2))
m6
1. 評価関数・目的関数¶
さて、最適を扱う分野であるオペレーションズリサーチの基本的な用語として解と評価関数（あるいは目的関数）があります。解とは、扱っている問題に対する答えです。今回は消防署の位置を答える問題ですので、消防署の位置が解に相当します。なお、解の良さはここでは議論していないことに注意してください。良い消防署の位置だろうが、悪い消防署の位置だろうが、解は解です。それで、その解の良さを測定するモノサシのことを、評価関数と呼びます。今回は評価関数として、市民一人当たりの消防署との平均距離を採用することにします。
- 解: 消防署5箇所の設置位置（経度・緯度5セット）
- 評価関数: 市民一人当たりの消防署との平均距離
評価関数が高い方がいいか低い方がいいかについては、扱う問題によって違います。今回の例で言えば当然低い方が良いことになります。が、「国民の幸せ」みたいなものを評価関数にした場合は、高い方が良いことになります。これを間違えないように注意しましょう。これを適当にやると、最良を目指していたつもりが、最悪を得ていた、ということになりかねません。
このために、評価関数をきちんと設計しなければなりません。この詳細は、以下の原稿にまとめられているので、それをご覧ください。
大前ほか, 遺伝的アルゴリズムによる千葉県習志野市の公共施設配置に対する近似解と施設削減, 第36回ファジィシステムシンポジウム, 2020.
普通にググればダウンロードできると思います。ない場合は、大前まで連絡のこと。
これを読むと、各座標と、その座標の人口と、消防署の座標、このデータが必要になることがわかります。これは、下記の変数として扱っています。参考までに、アクセスとそのデータ構造を表示してみましょう。
重要な変数：
- Place[i]: メッシュ地点iの座標（経度・緯度）
- Poplation[i]: メッシュ地点iの人口（人）
- FirePoint[i]: 消防署iの座標
NumMesh = len(Place)
NumFirePoint = len(FirePoint)
print(" +++ 基礎情報 +++")
print("メッシュ地点の総数: ", NumMesh)
print("消防署のの総数: ", NumFirePoint)
print("")
i=0
j=2
print("メッシュ地点", i, "の経度・緯度: ", Place[i])
print("メッシュ地点", i, "の人口: ", Population[i], " 人")
print("消防署", j, "の経度・緯度: ", FirePoint[j])
こんな感じで、メッシュ地点の座標と人口、現在設置されている消防署の座標にアクセスできます。
2. 遺伝的アルゴリズム¶
さて、次は良い消防署の配置を求めるアルゴリズムについて説明します。習志野市中に5箇所、全パターンに消防署を設置して、評価関数を求め、最も良い設置位置を発見することで、最適な位置を知ることが可能です。しかし、この方法では時間がいくらあっても計算が終わりません。
例えば、設置できるところが100箇所あったとして、その中から5箇所を選択して消防署を設置する組み合わせは、$_{100}C_{5} = 75,287,520$パターンもあります。 如何にコンピュータでも莫大な時間がかかってしまいますし、そもそも習志野市はすごく広いので100箇所どころではありません。
そのため、ある程度効率的に、最適とは言えないがより良い場所を探さなくてはなりません。この方法の一つに、遺伝的アルゴリズム（Genetic Algorithm; GA）があります。
GAの大雑把な流れは、まず、10〜100パターンくらい無作為に消防署5箇所の施設配置を決めて、その評価値を得ます。この中でいい配置を2箇所得ます。その2箇所の配置を色々と混ぜてみて、新しい配置場所の候補を10〜100パターン求めます。いい配置2箇所から多数の配置パターンを作ったので、たぶん、それらもいい配置です。...という処理を何度も繰り返して、より良い配置を探そうとします。ただしこの処理だと、習志野市内のごく一部しかみなくなってしまう場合があるので、時折、でたらめな配置パターンを混ぜ込んであげます。
なんとなくわかったところで、より詳細な説明を示します。
GAで登場する用語:
個体: 一つの解を指します。すなわち、5箇所の消防署の配置の1パターンの組み合わせを意味します。1箇所あたり、経度と緯度がありますので、1個体あたり$5 \times 2 = 10$個の数字で構成されます。現状存在する個体を現世代個体、それにより生成された次の個体を、次世代個体と呼びます。
遺伝子: 個体が有している一つ一つの情報を意味します。今回であれば1個体が10個の数字で構成されています。10個の数字の一つ一つを遺伝子と呼びます。
適応度（生存能力）: 個体が有している、評価関数の値を表します。生存能力が高い個体ほど、次世代に子孫を残せるという構造を持つのが生物ですが、この生存能力を評価関数と対応づけているわけです。
ランク戦略: 複数ある個体から、ある基準に基づいて2つの個体を取り出す操作を、選択と呼びます。この中で、各々の個体の評価値が高いほど、選択される確率が高いような選択方法を、ランク戦略と呼びます。例えば、個体A1, A2, A3, A4, A5がいたときに、その評価値を高い順に並べると、個体A2, A1, A4, A5, A3となったとします。この中から2つを選択するときに、A2を40%, A1を30%, A4を20%, A5を10%, A3を0%の確率で選択する場合、それはランク戦略を採用していることになります。なお、同じ個体は選択されないことにするのが普通です。他にも色々な選択方法があります。
一様交叉: 選択により選ばれた2つの個体において、それらの各遺伝子を一定確率で入れ替え、新たな個体を2セット生成する操作を一様交叉と呼びます。確率は50%が多いようです。個体A=abc、個体B=edfとあった場合に、確率サイコロを3回ふって、1桁目は入れ替える、2桁目を入れ替えない、3桁目を入れ替えない、そうなった場合は、個体C=ebc、個体D=adfという、新たな2つの個体が生成されます。
突然変異: ランク戦略により選択された1つ目の個体において、その遺伝子の1箇所に乱数で置き換えます（他にも色々な突然変異があります）。ただし、経度や緯度をランダムに生成すると、アメリカの変な街が選ばれるとか、かなり不都合があります。今回であれば習志野市内の経度緯度の最大最小を取得しておいて、その範囲内の乱数を生成することが必要でしょう。
遺伝的アルゴリズムにより最適施設の近似解を導出する手続き:
配置する施設数n、生成する個体数m、最大更新回数rを決定する。
メッシュ地点の経度・緯度の最大値と最小値を取得する。
その範囲内で、m個の個体を生成する。
生成されたm個の個体の評価関数に通し、評価値を得る。
ランク戦略により、2つの個体を取得する。
高確率で一様交叉、低確率で突然変異を実施し、個体を獲得する（99%, 1%など）。
生成される個体の数がmになるまで、処理4を繰り返す。
一様交叉では2個の個体が生成されるので、個体数がmを超える場合がある。この場合、個体数がmになるように個体を削除する。
最大更新回数に至っていないならば、処理3に進む。
最大更新回数に至ったならば、処理を終了する。
得られた個体の中で、評価値が最も良い解を「近似解」とする。
これが一連の手続きですが、以下のコードではより良い解を得るために、処理3に「m個の個体の評価値の標準偏差を算出し、標準偏差がほぼ0であれば、突然変異確率を上昇させる」という処理、大変異を入れています。全個体の評価値の標準偏差が0というのは、全個体がまったく同じ評価値を得ているので、すべてまったく個体であると考えられますから、交叉をしても無意味になります。よって、突然変異を行わないとより良い解を目指すことができなくなります。
選択や交叉は、ランク戦略や一様交叉以外にも他にも色々な種類があるので、興味がある人は調べてみましょう。
関数1. 評価関数¶
- 入力: メッシュ地点の人口と座標のリスト、1パターンの消防署配置（1個体のみ）
- 出力: その個体の評価値
入力が1個体のみですから、施設数を5としたとき、
FirePoint = [ [a, b], # 施設1の経度・緯度 [c, d], # 施設2の... [e, f], [g, h], [i, j] ]
こんな感じの、10個の数字を入れないと動かないので注意です。
# メッシュ地点の人口と座標のリスト、1パターンの消防署配置（1個体のみ）を入力
# その評価値を1つ返す関数
def Evaluater(Population, Place, FirePoint):
# +++ 人口メッシュ地点から、対応する消防署を一箇所のみ発見 +++ #
MinPoint = []
for i in range(len(Place)):
# メッシュ1地点と、消防署全地点の距離を求める
temp_list = []
for j in range(len(FirePoint)):
temp_list.append(distance.distance(Place[i], FirePoint[j]).km)
# 距離が最小の消防署のIDを取得
MinPoint.append(np.argmin(temp_list))
# +++ 評価関数の計算 +++ #
P = np.sum(Population) # 習志野市の総人口
PCX_each_i = []
for i in range(len(Place)):
# メッシュ地点iとそれをカバーする消防署との距離を算出
# x_ij=1の項のみを考慮して距離を算出しているので、sum c_ij * x_ij と等価
CX_i = distance.distance(Place[i], FirePoint[MinPoint[i]]).km
# CXに対し、人口をかける。総人口で割る
PCX_i = Population[i] * CX_i / P
PCX_each_i.append(PCX_i)
#print(PCX_each_i[i])
Eval = np.sum(PCX_each_i)
return Eval
# 使い方
Eval_sample = Evaluater(Population, Place, FirePoint)
print("習志野市の消防署に対する1人当たりの距離は", np.round(Eval_sample, 3), "kmです。")
関数2. 初期個体の生成¶
- 入力: 個体数、施設数、500mメッシュ座標
- 出力: 指定した個体数分の、個体ら
# 初期解候補の数、施設の数、500mメッシュ座標
def InitGenerator(gene_n, fire_n, Place):
import random
# 乱数生成範囲を取得 (人口メッシュの経度・緯度の最小最大)
point = np.zeros([len(Place), 2]) #経度・緯度2つぶん
for i in range(len(point[:, 0])):
for j in range(2):
point[i, j] = Place[i][j]
X_max = np.max(point[:,0])
X_min = np.min(point[:,0])
Y_max = np.max(point[:,1])
Y_min = np.min(point[:,1])
# 初期個体の生成
import random
Out_gene = []
for i in range(0, gene_n):
Each_Out_gene = []
for j in range(0, fire_n):
a = random.uniform(X_min, X_max)
b = random.uniform(Y_min, Y_max)
Each_Out_gene.append([a, b])
Out_gene.append(Each_Out_gene)
return Out_gene
# Out_gene[0]で0番目の解候補
# Out_gene[0][1]で0番目の解候補の、1番目の消防施設の経度緯度
# 使い方
InitGene = InitGenerator(3, 6, Place) # 生成する解候補の数, 施設数
print("0番目の解候補における、1番目の消防署の場所: ")
print(InitGene[0][1])
関数3. ランク戦略、一様交叉、突然変異¶
- 入力: m個の現世代の個体、事前個体の評価値、生成する個体数、メッシュ地点座標、突然変異確率、出力可視化（0ならしない、1ならする）
- 出力: m個の次世代個体
def MakeNextGene_Num(PreGeneList, EvalList, gene_n, Place, mutat_rate, progress):
import random
import copy
# 初期設定
th = mutat_rate # 突然変異確率（0〜1）
PostGene = [] # 次世代個体のリスト
# 個体生成スタート
while(len(PostGene) <= gene_n):
mutat_prob = random.uniform(0, 1) # 0-1サイコロを振る
# +++ 現世代個体を格納するゼロ行列を定義 +++ #
l_a = len(PreGeneList[0]) # 施設数
l_b = gene_n # 解候補の数
PreGeneNum_keido = np.zeros([l_b, l_a]) # 現世代個体（経度）
PreGeneNum_ido = np.zeros([l_b, l_a]) # 現世代個体（緯度）
PreGeneNum = np.zeros([l_b, l_a*2]) # 現世代個体（経度・緯度）
# +++ 現世代個体をnumpyへ代入 +++ #
# -> PreGeneNum[i, :]がi番目の解候補（施設数＊経度・緯度分の列サイズ）
for i in range(l_b):
for j in range(l_a):
PreGeneNum_keido[i, j] = PreGeneList[i][j][0]
PreGeneNum_ido[i, j] = PreGeneList[i][j][1]
for i in range(l_b):
k=0
for j in range(0, l_a*2, 2):
PreGeneNum[i, j] = PreGeneNum_keido[i, k]
PreGeneNum[i, j+1] = PreGeneNum_ido[i, k]
k=k+1
# +++ ランク選択 +++
if progress == 1:
print("+++ ランク選択・一様交叉 +++")
# -> 評価値が高い上位4つの個体を抽出
TempEvalList = copy.copy(EvalList) # id変更
r1 = np.argmin(EvalList)
TempEvalList[r1] = 9999 # 一番良い解の評価値を一旦下げる。
r2 = np.argmin(TempEvalList)
TempEvalList[r2] = 9999 # 一番良い解の評価値を一旦下げる。
r3 = np.argmin(TempEvalList)
TempEvalList[r3] = 9999 # 一番良い解の評価値を一旦下げる。
r4 = np.argmin(TempEvalList)
if progress == 1:
print("-> 上位4遺伝子: ", r1, r2, r3, r4)
# 確率による2遺伝子の抽出
prob_list = [0.4, 0.3, 0.2, 0.1] # 第1〜4位の個体選択確率（合計1）
same_flag=1
while(same_flag==1):
# 第1個体の抽出
Select_prob = random.uniform(0, 1) # 0-1サイコロを振る
if prob_list[0] > Select_prob: # 0~40%
s1 = r1
elif prob_list[0]+prob_list[1] > Select_prob: # 40~70%
s1 = r2
elif prob_list[0]+prob_list[1]+prob_list[2] > Select_prob: # 70~90%
s1 = r3
else: # 90~100%
s1 = r4
# 第2個体の抽出
Select_prob = random.uniform(0, 1) # 0-1サイコロを振る
if prob_list[0] > Select_prob: # 0~40%
s2 = r1
elif prob_list[0]+prob_list[1] > Select_prob: # 40~70%
s2 = r2
elif prob_list[0]+prob_list[1]+prob_list[2] > Select_prob: # 70~90%
s2 = r3
else: # 90~100%
s2 = r4
# 選択された個体が同じなら、再度選定
if progress == 1:
print("-> 仮選択個体: ", s1, s2)
if s1==s2:
same_flag=1
if progress == 1:
print("-> 選択遺伝子が一致したため、再度施行します。")
else:
same_flag=0
# 選択された2つの個体
SelectGene1 = PreGeneNum[s1, :] #個体1
SelectGene2 = PreGeneNum[s2, :] #個体2
if progress == 1:
print("-> 最終選択個体: ", s1, s2)
print("-> 遺伝子の一致率: ", np.round(np.sum(SelectGene1==SelectGene2)/len(SelectGene1), 2)*100, "%")
# 交叉 or 突然変異
if mutat_prob >= th:
# +++ 一様交叉 +++
# -> 上位2個体の各遺伝子を、50%の確率で入れ替えたものを次世代個体とする。
NewGene1 = copy.copy(SelectGene1) # id変更
NewGene2 = copy.copy(SelectGene2) # id変更
change_id = []
for i in range(0, l_a*2):
crossover_prob = random.uniform(0, 1) # 0-1サイコロを振る
if crossover_prob >= 0.5: # 50%ならRank1, 2のその値を入れ替え
NewGene1[i] = SelectGene2[i]
NewGene2[i] = SelectGene1[i]
change_id.append(i)
# 次世代個体の仮登録
PostGeneNum1=[]
PostGeneNum2=[]
for i in range(0, l_a*2, 2):
PostGeneNum1.append([NewGene1[i], NewGene1[i+1]])
PostGeneNum2.append([NewGene2[i], NewGene2[i+1]])
if progress == 1:
print("-> 一様交叉による入れ替えデータ: ", change_id)
# 次世代個体の本登録
PostGene.append(PostGeneNum1)
PostGene.append(PostGeneNum2)
if progress == 1:
print("-> 一様交叉により2個体が生成されました。")
print("")
# 突然変異
else:
if progress == 1:
print("+++ 突然変異 +++")
# 乱数生成範囲を取得 (人口メッシュの経度・緯度の最小最大)
point = np.zeros([len(Place), 2]) #経度・緯度2つ分
for i in range(len(point[:, 0])):
for j in range(2):
point[i, j] = Place[i][j]
X_max = np.max(point[:,0])
X_min = np.min(point[:,0])
Y_max = np.max(point[:,1])
Y_min = np.min(point[:,1])
# ランク戦略により選択された個体を取得
NewGene1 = copy.copy(SelectGene1) # id変更
insertID = random.randint(0, len(NewGene1)-1) #変更地点取得
if insertID%2 == 0: # 偶数なら、経度用の乱数生成
NewGene1[insertID] = random.uniform(X_min, X_max)
else: # 奇数なら、緯度用の乱数生成
NewGene1[insertID] = random.uniform(Y_min, Y_max)
# 次世代個体の仮登録
PostGeneNum1=[]
for i in range(0, l_a*2, 2):
PostGeneNum1.append([NewGene1[i], NewGene1[i+1]])
# 次世代個体の本登録
PostGene.append(PostGeneNum1)
if progress == 1:
print("-> 突然変異により個体が生成されました。")
print("")
# もし個体数がgene_nと等しくなければ、最後尾を削除
while(gene_n != len(PostGene)):
if progress == 1:
print("余分な個体を削除しました。")
PostGene.pop(-1)
if progress == 1:
print("")
return PostGene
メイン処理¶
以下、説明した処理通りにアルゴリズムを走らせる部分です。 交叉・突然変異の実行処理を見たい場合は、progress = 1にしてください。
# 事前準備
import time
import random
gene_n = 20# 生成する個体数（解候補の数）
fire_n = 5 # 配置する施設の数
max_step = 100 # 最大反復回数
progress = 0 # 実行処理を可視化しない
mutat_rate_low = 0.01 # 初期突然変異確率
mutat_rate_high = 0.50 # 解が安定した場合の突然変異確率
random.seed(19) # 乱数シード（再現可能な結果にしたい場合）
Gene = [] # 個体保存用リスト
Eval = [] # 評価値保存用リスト
# +++ 初期個体の生成 +++
Gene.append(InitGenerator(gene_n, fire_n, Place)) # 生成する解候補の数, 施設数
# 遺伝的アルゴリズム開始
for t in range(0, max_step):
start_time = time.time()
# 個体の評価値を取得
EachEval = []
for i in range(gene_n):
EachEval.append(Evaluater(Population, Place, Gene[t][i]))
#print(EachEval[i])
Eval.append(EachEval)
# 評価値のばらつきが極端に少ない場合は、突然変異確率を上昇させる。
if np.std(EachEval) < 0.001:
mutat_rate = mutat_rate_high
else:
mutat_rate = mutat_rate_low
# 次世代個体の生成
Gene.append(MakeNextGene_Num(Gene[t], EachEval, gene_n, Place, mutat_rate, progress))
end_time = time.time()
# 次世代遺伝子の登録
if t%int(max_step/5) == 0: # 指定回数の1/5進むたびに処理結果を表示
print("########## ", t, "世代 ##########")
print("")
print(" ********** 評価値の算出 **********")
print("-> 最小値: ", np.round(np.min(EachEval), 3))
print("-> 平均値: ", np.round(np.mean(EachEval), 3))
print("-> 標準偏差: ", np.round(np.std(EachEval), 3))
print("1処理に要する計算時間: ", np.round(end_time - start_time, 2), " [sec]")
print("")
print("終了！")
結果が終わったので、横軸を世代更新、縦軸を評価値として、学習状況を可視化してみます。
import matplotlib.pyplot as plt
PlotEval_mean = []
PlotEval_min = []
for i in range(max_step):
PlotEval_mean.append(np.mean(Eval[i])) # 第i世代個体の平均値
PlotEval_min.append(np.min(Eval[i])) # 第i世代個体の最小値
#plt.plot(PlotEval_mean, label="mean")
plt.plot(PlotEval_min, label="min")
plt.xlabel("Iteration")
plt.ylabel("Evaluation")
plt.legend()
世代更新されるたびに、評価値が減少していくことが確認できます。評価値は1人当たりの消防署との平均距離でしたから、市民と消防署の距離がどんどん近くなっているわけですね。うまくいかない場合は、乱数のシードを変えるなど、調整が必要です。
ところで、グラフをよくみてみると、10世代〜25世代や50〜85世代付近、ほどんと改善が見られないのに、それ以降は改善していることが見て取れます。これが突然変異の効果です。突然変異はランダムに優良な個体の1遺伝子を乱数に置き換えるので、悪い結果になってしまうことがほとんどですが、偶然、元の個体よりも良くなることがあります。これが、改善という結果を導くわけです。
良い解が得られたところで、それを地図上にプロットしてみます。
# 最終ステップの評価関数値Eval[-1]のなかで、最小のIDを取得
BestID = np.argmin(Eval[-1])
# 最終ステップの最適施設の座標を取得
BestFirePoint = Gene[-2][BestID]
# メッシュ地点に対する、利用する消防署IDを取得
MinPoint = DecisionPoint(BestFirePoint, Place)
# 習志野市を中心とするマップの作成
m7 = folium.Map(location=cent, zoom_start=13)
# 習志野市の塗りつぶし
folium.Choropleth(geo_data=geojson["features"][GetID],
fill_color='YlGn', # 塗りつぶしの色
fill_opacity=0.1, # 塗りつぶしの濃さ
line_opacity=0.9, # 境界の色の濃さ
).add_to(m7)
# 500mメッシュ地点のマーカを設定
MeanPop = np.mean(Population)
for i in range(0, len(Place)):
if Population[i] > MeanPop: # 平均オーバーなら、青
html_text = str(int(Population[i]))
folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_blue(html_text))).add_to(m7)
else:
html_text = str(int(Population[i]))
folium.Marker(location = Place[i], icon=DivIcon(html=MakeBox_red(html_text))).add_to(m7)
# 各消防署のカバー人口を算出
CoverPop = []
for i in range(len(BestFirePoint)):
s = 0
for j in range(len(Population)):
if MinPoint[j] == i: # そのメッシュ地点に対応している消防署のみをカウント
s = s + Population[j]
CoverPop.append(s) # 消防署iが対応する総人数をappend
# 消防署をプロット
for i in range(len(FireName)):
folium.Marker(location = BestFirePoint[i],
popup = "<p>Fire station: "+ str(i)+"</p><p>Cover: " + str(int(CoverPop[i])) + " persons</p>").add_to(m7)
# 消防署と500mメッシュのラインを表示
for i in range(len(Place)):
folium.PolyLine(locations=[BestFirePoint[MinPoint[i]], Place[i]], weight=2, color='black').add_to(m7)
# 凡例はhtmlで入れるしかない模様
legend_html = '''
<div style="position: fixed;
top: 20px; left: 50px; width: 210px; height: 67px; Background: white;
border:2px solid grey; z-index:9999; font-size:14px; padding-left: 5px">
Blue box: High population area.<br>
Red box: Low population area.<br>
Marker: Fire station
</div>
'''
legend_html2 = '''
<div style="position: fixed;
top: 1px; left: 50px; width: 260px; height: 67px; Background: transparent;;
border:0px solid grey; z-index:9999; font-size:14px; padding-left: 5px">
Made by Yuto Omae at Nihon Univ.
</div>
'''
m7.get_root().html.add_child(folium.Element(legend_html))
m7.get_root().html.add_child(folium.Element(legend_html2))
m7
EvalBestPoint = Evaluater(Population, Place, BestFirePoint)
print("習志野市の消防署に対する1人当たりの距離は", np.round(EvalBestPoint, 3), "kmです。")
GAにより得られた配置の場合、消防署と習志野市民の1人当たりの平均距離は、683mだそうです。現在の配置との距離を比較すると、以下のようになります。
- 実際の施設配置: 828 m
- GAによる施設配置: 683 m
150mくらい近くすることができました。市民全員の平均値なわけですから、かなり良い効果をもたらしたことになります。もちろん、土地の代金も違いますし、すでに建物もありますし、消防法上ここに建てちゃダメとかそういう難しいルールもあるでしょう。故に、そこに消防署を建てろといっても土台無理話ではあります。が、今回の事例紹介を厳密に理解して、自分で変更を加えられるくらいに十分なプログラミング能力があれば、以下のように現実的に使えそうな問題にも対応することが可能です。
- ラーメン屋さんを出店するとき、競合店と最も距離が遠くなる、かつ、市民が多い、これらを満たす配置はどこか。
- 消防署を一つ減らさなければならないという要求があったときに、どこを減らすことが良いか。
- など...
都市計画など、大きな仕事をやる場合に必要な能力とも言えるかもしれません。
保存しておしまい。
m7.save("OptFire.html")
htmlファイルなので、もちろん、httpのサーバにおけば外部公開できます！！
本当の終わり。