蟻コロニー最適化: サラリーマンが土曜日の昼下がりにセールスマンと蟻で遊んでたら夜になった話

読書の秋、最適化の秋ですね。

今日はサラリーマンが休日にセールスマンと蟻とで遊んだ記録です。
ぜんぜん寂しくないです。最適な秋の過ごし方の一つです。

サラリーマンが、セールスマンと蟻で遊んだ https://t.co/pzGyKE84QB

— えんせき@ノースリーブタートルネック (@EnsekiTT) 2016年10月30日

赤い線はフェロモンで色の濃さはフェロモンの濃さです。 青い線は実際にある1匹のありが通ったところです。

蟻の戦略

http://science.newsln.jp/articles/2014052623200004.html
・蟻が適当に餌を探して歩き回る
・餌を見つけると道にフェロモンを残しながら巣に戻る
・餌を見つけられず疲れると巣に戻る
・次の蟻はフェロモンを目印にしつつ餌を探して歩き回る
ここでは一匹の蟻に注目しましたが、これをうじゃうじゃ居る蟻たちがそれぞれ行います。

巡回するセールスマンの職場環境

http://elib.zib.de/pub/Packages/mp-testdata/tsp/tsplib/tsplib.html
・セールスマンは各都市を1回ずつ訪問しなければいけない
・セールスマンは各都市間の移動時間を知っている
・セールスマンは可能な限り早く仕事を終わらせたい
2,3箇所なら問題ない職場環境でしょう、しかしこれが20箇所を超えてくると大変なことになります。

問題設定

セールスマン「困難💢」ってことで蟻の力を借りることに。
訪問場所にそれぞれ餌を置き、街中に蟻を放ち、それらを結ぶ最短経路を求めてもらうことにしました。
訪問される側からしたらたまったものではありませんね。

アルゴリズム

巡回セールスマン問題を蟻コロニー最適化で最適解っぽいのを探します。

ウィキペディアに書いてあったアルゴリズム

エージェント（アリ）とフェロモンの初期化
終了条件を満たすまで以下の処理を繰り返す。
各エージェントに対して、フェロモンとヒューリスティックな情報に基づいて確率的な解の選択を行う。
各エージェントが分泌するフェロモンを計算する。
フェロモン情報の更新
最も良い成績のエージェントの解を出力する。

結局なんのかんのやって実際に実装したアルゴリズム

都市の生成（訪問先の数と訪問先間の距離を決める）
イテレーションの数だけ以下の処理を繰り返す
　a. 各蟻を生成する（今いる訪問先と訪問先間のフェロモンの情報と訪問した先、してない先のリストを与える）
　b. 各蟻がすべてを訪問するまで以下の処理を繰り返す
　　ア. 現在の訪問先から次の訪問先候補の評価値を決める
　　イ. 各訪問先に蟻が行く確率を求める（蟻は確率を考慮しつつも割と適当に動き回る）
　　ウ. 確率に応じて移動する
　c. フェロモンを更新する
3．最も最短時間で全都市回れたエリート蟻が決まる

記号はWikipedia準拠で
https://ja.wikipedia.org/wiki/%E8%9F%BB%E3%82%B3%E3%83%AD%E3%83%8B%E3%83%BC%E6%9C%80%E9%81%A9%E5%8C%96
最終更新 2016年10月15日 (土) 19:18 （日時は個人設定で未設定ならばUTC）。
のもの

評価値

a_{{ij}}^{{k}}(t)={\frac {[\tau _{{ij}}(t)]^{{\alpha }}[\eta _{{ij}}]^{{\beta }}}{\sum _{{l\in \Omega }}[\tau _{{il}}(t)]^{{\alpha }}[\eta _{{il}}]^{{\beta }}}}

移動確率

p_{{im}}^{{k}}(t)={\frac {a_{{im}}^{{k}}(t)}{\sum _{{n\in \Omega }}a_{{in}}(t)}}

フェロモン分泌量

\Delta \tau _{{ij}}^{{k}}(t)={\begin{cases}Q/L_{k}(t),&{\mbox{if }}(i,j)\in T_{k}(t)\\0,&else\end{cases}}

フェロモンの更新（Wikipediaとここだけ変えた）

\tau _{{ij}}(t+1)=\rho \tau _{{ij}}(t)+(1-\rho)\sum _{{k=1}}^{{m}}\Delta \tau _{{ij}}^{{k}}(t)

※フェロモンの更新が収束したらフェロモンを初期化して局所解から抜けられる仕組みを追加

変数:

訪問先の数: 30
蟻の数: 20
蟻を放つ回数: 500
スタート地点: 0
α(フェロモンの優先度): 1
β(距離の優先度): 5
Q(総距離の影響を左右する定数): 100
ρ(フェロモンの蒸発率): 0.4
1-ρ(フェロモンの定着率)

環境:

可視化のためにMatplotlibを入れたPython3.5
実行するディレクトリにhistという蟻の軌跡を保存するためのディレクトリを用意しておく
（下のコードでは生成しない）

コード:

蟻（エージェント）のクラス

class Agent():
    def __init__(self, towns, roads, start, pheromone):
        # value
        self.current = start
        self.alpha = 1
        self.beta = 3
        self.Q = 100

        # list
        self.whole = copy.deepcopy(towns)
        self.notVisited = copy.deepcopy(towns)
        self.notVisited.remove(start)
        self.visited = [towns[start]]

        # dict
        self.roads = roads
        self.pheromone = copy.deepcopy(pheromone)
        self.way = {}
        for i in towns:
            for j in towns:
                if i is j:
                    continue
                self.way[(i,j)] = False

    def assessment(self, j):
        numerator = (self.pheromone[(self.current, j)]**self.alpha) * ((1/self.roads[(self.current, j)])**self.beta)
        denominator = sum([(self.pheromone[(self.current,l)]**self.alpha) * ((1/self.roads[(self.current, l)])**self.beta) for l in self.notVisited])
        assess = numerator / denominator
        return assess

    def probability(self):
        p = {}
        for m in self.notVisited:
            mAsses = self.assessment(m)
            sumAsses = sum([self.assessment(n) for n in self.notVisited])
            p[m] = mAsses / sumAsses
        return p

    def agentwalk(self):
        for i in self.whole:
            prob = self.probability()
            sprob = prob.items()
            #sprob = sorted(prob.items(), key=lambda x: x[0])
            choice = random.random()
            for i in sprob:
                choice = choice - i[1]
                if choice < 0:
                    nextT = i
                    break
            self.way[(self.current, nextT[0])] = True
            self.current = nextT[0]
            self.visited.append(nextT[0])
            self.notVisited.remove(nextT[0])
            if len(self.notVisited) == 0:
                return self.visited
        return self.visited

    def get_deltaphero(self):
        sumoflen = self.get_length()
        deltaPheromone = {}
        for i in self.pheromone:
            if self.way[i]:
                deltaPheromone[i] = self.Q/sumoflen
            else:
                deltaPheromone[i] = 0
        return deltaPheromone

    def get_phero(self,startpheromone, pheromones, ro):
        for i in startpheromone:
            startpheromone[i] = ro * startpheromone[i] + (1-ro) * sum([k[i] for k in pheromones])
        return startpheromone

    def get_length(self):
        sumoflen = 0
        for i in range(1,len(self.visited)):
            sumoflen += roads[(self.visited[i-1],self.visited[i])]
        return sumoflen

    def get_way(self):
return self.visited

可視化のコード

figureCount = 0
def anthistorySave(len, way, delt, pheromone, title=""):
    global figureCount
    phrange = max(pheromone.values()) - min(pheromone.values()) + 1
    for i in pheromone:
        pheromone[i] = pheromone[i]/phrange

    fig, ax = plt.subplots()
    plt.title(title + " length=" + str(len))
    plt.scatter([i[0] for i in positions], [i[1] for i in positions])
    for i, town in enumerate(way):
        plt.annotate(i, (positions[town][0]+1, positions[town][1]-3), color='r')
        plt.annotate(town, (positions[town][0]+1, positions[town][1]+1), color='g')
    for i in pheromone:
        plt.plot([positions[i[0]][0], positions[i[1]][0]],
                 [positions[i[0]][1], positions[i[1]][1]], 'r', alpha=pheromone[i])
        if delt[i] > 0:
            plt.plot([positions[i[0]][0], positions[i[1]][0]],
                     [positions[i[0]][1], positions[i[1]][1]], 'b')
    plt.grid()
    plt.savefig("hist/" + title + "anthistory_"+"{0:05d}".format(figureCount)+".png")
    plt.close('all')
figureCount += 1

実行するコード

if __name__ == "__main__":
    NumOfTown = 30
    NumOfAgent = 20
    NumOfSolve = 500
    ro = 0.4
    random.seed(2)

    towns = [i for i in range(NumOfTown)]
    positions = [(random.randint(1,100), random.randint(1,100)) for i in range(NumOfTown)]

    roads = {}
    pheromone = {}
    for i in towns:
        for j in towns:
            if i is j:
                continue
            roads[(i,j)] = sqrt((positions[i][0] - positions[j][0])**2 +
                            (positions[i][1] - positions[j][1])**2)
            pheromone[(i,j)] = random.random()

    minlength = None
    bestAgent = None
    lastPheno = 0
    for i in range(NumOfSolve):
        pheros = []
        topAgent = None
        for m in range(NumOfAgent):
            k = Agent(towns=towns, start=0, roads=roads, pheromone=pheromone)
            k.agentwalk()
            length = k.get_length()
            pheros.append(k.get_deltaphero())
            if (topAgent is None) or (topAgent.get_length() > length):
                topAgent = copy.deepcopy(k)
            if (minlength is None) or (minlength > length):
                minlength = length
                bestAgent = copy.deepcopy(k)
            anthistorySave(k.get_length(), k.get_way(), k.get_deltaphero(), k.pheromone)
        #anthistorySave(topAgent.get_length(), topAgent.get_way(), topAgent.get_deltaphero(), topAgent.pheromone)
        print("今" + str(i) + "番目のありんこたちが仕事をしました。")
        pheromone = k.get_phero(pheromone, pheros, ro)

        if round(sum(pheromone.values()),4) == round(lastPheno,4):
            pheromone = {}
            for i in towns:
                for j in towns:
                    if i is j:
                        continue
                    pheromone[(i,j)] = random.random()
        lastPheno = sum(pheromone.values())

    print(bestAgent.get_way())
anthistorySave(bestAgent.get_length(), bestAgent.get_way(), bestAgent.get_deltaphero(), bestAgent.pheromone)

今後の展望

次の休暇は焼きなまして遊ぼうかと思う。

クリエイティブ・コモンズ 表示-継承ライセンス