2014年12月6日土曜日

pythonで遺伝的アルゴリズムの練習2 --巡回セールスマン問題

前回に引き続き、遺伝的アルゴリズムの例題をもう一つ解いてみようと思います。
今回は巡回セールスマン問題を解きます。
巡回セールスマン問題とは、例えば5つの都市があったとして、すべて都市を1回づつ訪問しようと思っている場合、
どのような順番で都市を回れば最短ルートで行くことができるのかという問題です。
都市の数が5つぐらいなら総当りでなんとかなるのですが、30、50となってくると、
組み合わせが階乗になるのでとても総当りできません。
そこで、この問題を解くのに様々な方法があるのですが、その方法の1つとして遺伝的アルゴリズムが使われています。
遺伝的アルゴリズムでは、最適解を出すことは難しく、
最適解に限りなく近い解しか出すことができませんが、いい練習になると思います。



まず、遺伝子ですが、都市のリストが最初から与えられているとして、
都市に0から順番に番号を振ります。
そして、都市番号を遺伝子配列一つ一つにランダムに与えます。
この時一度使った数字は再登場させないようにします。
つまり、都市[A,B,C,D,E]とある場合、
A→0
B→1
.
.
E→4
と都市番号を割り当てます。
そして、遺伝子情報としては[0,1,2,3,4]が[A→B→C→D→E]となり、
[3,0,4,1,2]が[D→A→E→B→C]となります。
このようにして書いてみようと思ったのですが、
このままでは交叉や突然変異をした時に、同じ都市にもう一度訪問してしまう可能性があるので、
実際に遺伝子にコードする情報は、[ABCDE]の順に都市を回るなら、[0,0,0,0,0]
[DAEBC]の順なら、[3,0,2,0,0]のように書くことにします。
これは一度訪問した都市を除外した未訪問都市リストを作り、そのリストの中の何番目かを遺伝子にコードしています。
このような書き方を順序表現というらしいです。
また、最初の遺伝子のまま先に進めていく方法もあるようなのですが、こちらの方では交叉や突然変異の時に少し工夫が必要なようです。

また、評価方法は総距離のみで評価することにしました。

# -+- coding: utf-8 -*-
import random #ランダムモジュール
import math
import copy
import numpy as np
import pylab

def city_init():
    #都市座標
    citylist_x = [2,10,7,32,35,9,28,3,8,5,13,45,63,11,17,34,55,38,41,15,6,53,46,18,23,16,19,27,42,12]
    citylist_y = [36,22,47,41,32,24,15,2,25,28,48,5,11,16,27,8,46,35,58,14,26,21,49,34,1,7,19,43,29,39]
    return [citylist_x,citylist_y]

##評価関数
def eval_func(gean):
    x,y =city_init()
    x_cp = list(x)
    y_cp = list(y)
##    print x
    #a,b,c,d,e,f,のうち行った都市は消していく
    #c,e,a,f,d,bの順で回る場合は
    #[2,3,0,2,1,0]となる
    ##gean = [2, 3, 1, 1, 0]
    ##gean = [2,4,1,3,0]
    gean_decode = []
    gean_number = 0
    while gean_number<=len(gean)-1:
        num = x.index(x_cp[gean[gean_number]])
        gean_decode.append(num)
        del x_cp[gean[gean_number]]
        gean_number +=1
    #評価の基準は距離が短いことである。
    dist_stok = []
    gean_number = 0
    while gean_number<=len(gean)-1:
        dist = (x[gean_decode[gean_number+1]]-x[gean_decode[gean_number]])**2 + (y[gean_decode[gean_number+1]]-y[gean_decode[gean_number]])**2
        dist_stok.append(dist)
        gean_number +=1
    alldistance = sum(dist_stok)

    last_dist = (x[gean_decode[len(gean)-1]]-x[gean_decode[0]])**2 + (y[gean_decode[len(gean)-1]]-y[gean_decode[0]])**2
    alldistance = alldistance + last_dist
##    print alldistance
    return alldistance


def geneticoptimize(maxiter = 300000,maximize = False,popsize = 10000,popnum = 30,elite = 0.99,mutprob =0.3):
    """
    maxiter = 1, 繰り返し数
    maximize = True,    スコアを最大化
    popsize = 50,   個体数
    popnum = 10,    遺伝子数(長さ)
    elite = 0.2,    生き残る遺伝子の割合
    mutprob =0.2    突然変異のおこる確立
    """
    # 突然変異
    def mutate(vec):
        i = random.SystemRandom().randint(0,popnum-1)
        after_mut = random.randint(0,popnum-i-1)
        return vec[:i] + [after_mut]+vec[i+1:]

     # 1点交叉 非推奨
    def one_point_crossover(r1,r2):
        i = random.SystemRandom().randint(1,popnum-2)

        return random.SystemRandom().choice((r1[0:i] + r2[i:], r2[0:i] + r1[i:]))

    # 2点交叉
    def two_point_crossover(r1,r2):
        i, j = sorted(random.SystemRandom().sample(range(popnum),2))
        return random.SystemRandom().choice((r1[0:i] + r2[i:j] + r1[j:] , r2[0:i] + r1[i:j] + r2[j:]))

    # 一様交叉
    def uniform_crossover(r1, r2):
        q1 = copy.copy(r1)
        q2 = copy.copy(r2)
        for i in range(len(r1)):
            if random.SystemRandom().random() < 0.5:
                q1[i], q2[i] = q2[i], q1[i]

        return random.SystemRandom().choice([q1,q2])
    #遺伝子の初期化
    pop = []
    vec = []
    for i in range(popsize):
        selected_city = 1
        while selected_city <= popnum:
            vec_vec = random.SystemRandom().randint(0,popnum-selected_city)
            vec.append(vec_vec)
            selected_city +=1
        pop.append(list(vec))   #listはコピーを渡さないと変更元まで変更される
        del vec[:]


    #遺伝子の初期化

    # 交叉アルゴリズムの選択
    crossover = two_point_crossover
    #crossover = uniform_crossover

    #メインループ
    topelite = int(elite * popsize)
    score_stoc = [0]
    for i in range(maxiter):
        scores=[(eval_func(v),v) for v in pop]
        scores.sort()
        if maximize:
            scores.reverse()
        ranked = [v for (s,v) in scores]
        # 弱い遺伝子は淘汰される
##        pop = ranked[0:topelite]
        eliminate_cycle = 0

        while eliminate_cycle < popsize:
            if eliminate_cycle == 0:
                pop = ranked[:1]
            else:
                if random.SystemRandom().random() < elite**eliminate_cycle:
                    pop.append(ranked[eliminate_cycle])
            eliminate_cycle +=1

        while len(pop) < popsize:
            if random.SystemRandom().random() < mutprob:
                # 突然変異
                c = random.SystemRandom().randint(0,topelite)
                pop.append(mutate(ranked[c]))

            else:
                # 交叉
                c1 = random.SystemRandom().randint(0,topelite)
                c2 = random.SystemRandom().randint(0,topelite)
                pop.append(crossover(ranked[c1],ranked[c2]))
        # 暫定の値を出力
        #print(scores[0][0])
        #print(scores[0])
        score_stoc.append(scores[0][0])
    return scores[0]


def main():
    x,y =city_init()
    ans = geneticoptimize()
    ans_x = []
    ans_y = []
    print("Ans:",ans)
    for i in ans[1]:
        ans_x.append(x[i])
        ans_y.append(y[i])
        del x[i],y[i]
    print ans_x
##    ans_x.sort()
    print ans_x
    pylab.plot(ans_x,ans_y)
    pylab.show()



if __name__ == '__main__':
    main()

となります。
また、淘汰の条件も変更しました。
前回は、ソートした順に何番目までが生存としていたのですが、
これだと有料遺伝子を持っていたとしても他の因子により総距離が伸びてしまうと問答無用で淘汰されていたのですが、
今回は1位は100%2位は0.97%3位は0.97^2%4位は0.97^3というようにしました。
これで下位のものでも生存する確立が上昇します。


参考

http://aidiary.hatenablog.com/entry/20110109/1294581648



0 件のコメント:

コメントを投稿