手法としてはいろいろあると思うのですが，今回は遺伝的アルゴリズム(GA)と蟻コロニー最適化(ACO)の2つを実装・比較しました．
1.巡回セールスマン問題
よく出される例としてはゴミ収集車の経路でしょうか．ゴミ収集所から出発して，100軒の家を最短で全て回ってゴミ収集所にまた戻ってくるまでの経路を求める話です．
まぁ僕がここで云々説明せずともネット上には沢山関連するサイトが有りますからそちらを参照してください．とりあえずWikipediaと検索したらトップに出てきたサイトのURLを貼っときます．
真ん中のサイトは手法について非常に詳しく解説されている上，Pythonで実装されたソースコードがありますのでそちらもごらんになると良いかと思います．
一番下のサイトは巡回セールスマンについての説明が詳しく，フリーソフトの紹介もされています．
この問題はまだ多項式時間で最適解を解くことができるアルゴリズムが無く，うちの大学でも研究してる研究室がたしかあった気がします．
2.蟻コロニー最適化による実装
蟻コロニー最適化に関する数式は今回全てWikipediaに載っていたものを利用しました（蟻コロニー最適化 - Wikipedia）．
一応『心はプログラムできるか』という本を読んだ時に出てきてまして非常にわかりやすい解説がなされていましたので読んでみると良いかと思います．
以下ソースコード
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #coding:utf-8 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import pandas as pd | |
| from scipy.spatial import distance as dis | |
| """ | |
| 参考URL | |
| [1] 蟻コロニー最適化 - Wikipedia https://ja.wikipedia.org/wiki/蟻コロニー最適化 | |
| [2] 任意の確率密度分布に従う乱数の生成（von Neumannの棄却法） | Pacocat's Life http://pacocat.com/?p=596 | |
| """ | |
| class TSP: | |
| def __init__(self,path=None,alpha = 1.0,beta = 1.0,Q = 1.0,vanish_ratio = 0.95): | |
| """ 初期化を行う関数 """ | |
| self.alpha = alpha # フェロモンの優先度 | |
| self.beta = beta # ヒューリスティック情報(距離)の優先度 | |
| self.Q = Q # フェロモン変化量の係数 | |
| self.vanish_ratio = vanish_ratio # 蒸発率 | |
| if path is not None: | |
| self.set_loc(np.array(pd.read_csv(path))) | |
| def set_loc(self,locations): | |
| """ 位置座標を設定する関数 """ | |
| self.loc = locations # x,y座標 | |
| self.n_data = len(self.loc) # データ数 | |
| self.dist = dis.squareform(dis.pdist(self.loc)) # 距離の表を作成 | |
| self.weight = np.random.random_sample((self.n_data,self.n_data)) # フェロモンの量 | |
| self.result = np.arange(self.n_data) # もっともよかった順序を保存する | |
| def cost(self,order): | |
| """ 指定された順序のコスト計算関数 """ | |
| n_order = len(order) | |
| return np.sum( [ self.dist[order[i],order[(i+1)%n_order]] for i in np.arange(n_order) ] ) | |
| def plot(self,order=None): | |
| """ 指定された順序でプロットする関数 """ | |
| if order is None: | |
| plt.plot(self.loc[:,0],self.loc[:,1]) | |
| else: | |
| plt.plot(self.loc[order,0],self.loc[order,1]) | |
| plt.show() | |
| def solve(self,n_agent=1000): | |
| """ 巡回セールスマン問題を蟻コロニー最適化で解く """ | |
| order = np.zeros(self.n_data,np.int) # 巡回経路 | |
| delta = np.zeros((self.n_data,self.n_data)) #フェロモン変化量 | |
| for k in range(n_agent): | |
| city = np.arange(self.n_data) | |
| now_city = np.random.randint(self.n_data) # 現在居る都市番号 | |
| city = city[ city != now_city ] | |
| order[0] = now_city | |
| for j in range(1,self.n_data): | |
| upper = np.power(self.weight[now_city,city],self.alpha)*np.power(self.dist[now_city,city],-self.beta) | |
| evaluation = upper / np.sum(upper) # 評価関数 | |
| percentage = evaluation / np.sum(evaluation) # 移動確率 | |
| index = self.random_index(percentage) # 移動先の要素番号取得 | |
| # 状態の更新 | |
| now_city = city[ index ] | |
| city = city[ city != now_city ] | |
| order[j] = now_city | |
| L = self.cost(order) # 経路のコストを計算 | |
| # フェロモンの変化量を計算 | |
| delta[:,:] = 0.0 | |
| c = self.Q / L | |
| for j in range(self.n_data-1): | |
| delta[order[j],order[j+1]] = c | |
| delta[order[j+1],order[j]] = c | |
| # フェロモン更新 | |
| self.weight *= self.vanish_ratio | |
| self.weight += delta | |
| # 今までで最も良ければ結果を更新 | |
| if self.cost(self.result) > L: | |
| self.result = order.copy() | |
| # デバッグ用 | |
| print("Agent ... %d , Cost ... %lf" % (k,self.cost(self.result))) | |
| return self.result | |
| def random_index(self,percentage): | |
| """ 任意の確率分布に従って乱数を生成する関数 """ | |
| n_percentage = len(percentage) | |
| while True: | |
| index = np.random.randint(n_percentage) | |
| y = np.random.random() | |
| if y < percentage[index]: | |
| return index | |
| if __name__=="__main__": | |
| #tsp = TSP(path="input.csv") | |
| tsp = TSP() | |
| tsp.set_loc(np.random.random_sample((30,2))) | |
| tsp.plot() # 計算前 | |
| tsp.solve(n_agent=1000) # 1000匹の蟻を歩かせる | |
| tsp.plot(tsp.result) # 計算後 | |
3.遺伝的アルゴリズムによる実装
ぶっちゃけ蟻コロニー最適化よりも精度低いです．30個点でさえも蟻コロニー最適化に負けています．1000個ぐらいの点になるともう最適解求めるなんて夢のまた夢でした．原因としては恐らく突然変異がランダムすぎて，ほとんど改悪になるばかりという結果でしょうか．遺伝子は巡回経路の順番を表した変数(5つの点なら[2,4,3,1,0]という配列)です．突然変異で中身を入れ替えてます．
※追記：先日図書館にてほかの人の実装例を見ましたがこれとは異なる手法を用いていました．それを踏まえるとこのプログラムは一般的な実装例とは異なると思われますので注意してください．
ソースコード
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| #coding:utf-8 | |
| import numpy as np | |
| import matplotlib.pyplot as plt | |
| import pandas as pd | |
| from scipy.spatial import distance as dis | |
| class TSP: | |
| def __init__(self,path=None,n_gene = 256,n_parent = 10,change_ratio = 0.1): | |
| """ 初期化を行う関数 """ | |
| self.n_gene = n_gene # 一世代の遺伝子の個数 | |
| self.n_parent = 10 # 親として残す個体数 | |
| self.change_ratio = change_ratio # 突然変異で変化させる場所の数 | |
| if path is not None: | |
| self.set_loc(np.array(pd.read_csv(path))) | |
| def init_genes(self,): | |
| """ 遺伝子をランダムに初期化 """ | |
| self.genes = np.zeros((self.n_gene,self.n_data),np.int) | |
| order = np.arange(self.n_data) | |
| for i in range(self.n_gene): | |
| np.random.shuffle(order) | |
| self.genes[i] = order.copy() | |
| self.sort() | |
| def set_loc(self,locations): | |
| """ 位置座標を設定する関数 """ | |
| self.loc = locations # x,y座標 | |
| self.n_data = len(self.loc) # データ数 | |
| self.dist = dis.squareform(dis.pdist(self.loc)) # 距離の表を作成 | |
| self.init_genes() # 遺伝子を初期化 | |
| def cost(self,order): | |
| """ 指定された順序のコスト計算関数 """ | |
| return np.sum( [ self.dist[order[i],order[(i+1)%self.n_data]] for i in np.arange(self.n_data) ] ) | |
| def plot(self,order=None): | |
| """ 指定された順序でプロットする関数 """ | |
| if order is None: | |
| plt.plot(self.loc[:,0],self.loc[:,1]) | |
| else: | |
| plt.plot(self.loc[order,0],self.loc[order,1]) | |
| plt.show() | |
| def solve(self,n_step=1000): | |
| """ 遺伝的アルゴリズムで解く関数 """ | |
| for i in range(n_step): | |
| print("Generation ... %d, Cost ... %lf" % (i,self.cost(self.genes[0]))) | |
| self.step() | |
| self.result = self.genes[0] | |
| return self.result | |
| def step(self): | |
| """ 遺伝子を一世代分進化させる関数 """ | |
| # 突然変異 | |
| for i in range(self.n_parent,self.n_gene): | |
| self.genes[i] = self.mutation( np.random.randint(self.n_parent) ) | |
| self.sort() # ソートする | |
| def sort(self): | |
| """ 遺伝子を昇順にする関数 """ | |
| # コストを計算し，ソート | |
| gene_cost = np.array([self.cost(i) for i in self.genes]) | |
| self.genes = self.genes[ np.argsort(gene_cost) ] | |
| def mutation(self,index): | |
| """ 突然変異を起こした遺伝子を返す関数 """ | |
| n_change = int(self.change_ratio * self.n_data) | |
| gene = self.genes[index].copy() | |
| for i in range(n_change): | |
| # n_changeの個数だけ値を入れ替える | |
| left = np.random.randint(self.n_data) | |
| right = np.random.randint(self.n_data) | |
| temp = gene[left] | |
| gene[left] = gene[right] | |
| gene[right] = temp | |
| return gene | |
| if __name__=="__main__": | |
| #tsp = TSP(path="input.csv") | |
| tsp = TSP() | |
| tsp.set_loc(np.random.random_sample((30,2))) | |
| tsp.plot() | |
| tsp.solve() | |
| tsp.plot(tsp.result) |
4.評価
一応同じデータを用いてやった結果です．30個の[0,1)のランダムな二次元の点について行いました．パラメータはデフォルトです．まずは遺伝的アルゴリズム．
直線が交差してますね．コストは6.056765でした．続いて蟻コロニー最適化．
とりあえず今回の実装では蟻コロニー最適化のほうが早かったですが，遺伝的アルゴリズムも突然変異をもう少し工夫してやればもっと良い結果が出るような気がします．いつか焼きなまし法とかやってみたいです．
